WARCL polar map to energy-momentum mapping

// $BACKGROUND$

Number Ne = 670/1 // Number of energy bins
Number Nm = 320*4 // Number of pixels in momentum map
Number Nphi = 180 // Number of phi slices over 180 deg
Number MapRotation = 0

Image raw := GetFrontImage()
Number Nx = raw.ImageGetDimensionSize(0)
Number Ny = raw.ImageGetDimensionSize(1)
Number Nz = raw.ImageGetDimensionSize(2)
Number Dispersion_nm = raw.ImageGetDimensionScale(2)
Number WavelengthMin_nm = raw.ImageGetDimensionOrigin(2)
Number WavelengthMax_nm = WavelengthMin_nm + Dispersion_nm*(Nz -1)
Number WavelengthCenter_nm = WavelengthMin_nm + Dispersion_nm*(Nz-1)/2
Number EnergyMin_eV = 1239.8/(WavelengthMax_nm)
Number EnergyMax_eV = 1239.8/(WavelengthMin_nm)
Number Kmax_um = sin(90*pi()/180)*2*pi()/(WavelengthMin_nm/1000)
Number AngleSampling_deg = raw.ImageGetDimensionScale(0)
Number AngleMax_deg = AngleSampling_deg*Nx/2

Result("\n" + Dispersion_nm)
Result("\n" + WavelengthMin_nm)
Result("\n" + WavelengthCenter_nm)
Result("\n" + EnergyMin_eV)
Result("\n" + EnergyMax_eV)
Result("\n" + AngleSampling_deg)
Result("\n" + AngleMax_deg)

Image Energy_eV := RealImage("Energy", 4, Nx, Ne)
Energy_eV = EnergyMin_eV + (EnergyMax_eV - EnergyMin_eV)*irow/Ne

Image Row := (Nz/(WavelengthMax_nm - WavelengthMin_nm))*(1239.8/Energy_eV - WavelengthMin_nm)
//Row.showImage()

Image RawSlice := RealImage("Temp raw slice", 4, Nx, Nz)
Image EnergyARCL := RealImage("Energy sampled ARCL", 4, Nx, Ny, Ne)

EnergyARCL.ImageSetDimensionScale(2, (EnergyMax_eV - EnergyMin_eV)/Ne)
EnergyARCL.ImageSetDimensionOrigin(2, EnergyMin_eV)
EnergyARCL.ImageSetDimensionUnitString(2, "eV")
EnergyARCL.ImageSetDimensionScale(0, raw.ImageGetDimensionScale(0))
EnergyARCL.ImageSetDimensionOrigin(0, raw.ImageGetDimensionOrigin(0))
EnergyARCL.ImageSetDimensionUnitString(0, raw.ImageGetDimensionUnitString(0))
EnergyARCL.ImageSetDimensionScale(1, raw.ImageGetDimensionScale(1))
EnergyARCL.ImageSetDimensionOrigin(1, raw.ImageGetDimensionOrigin(1))
EnergyARCL.ImageSetDimensionUnitString(1, raw.ImageGetDimensionUnitString(1))

//EnergyARCL.showImage()

for (number i = 0; i < Ny; i++)
{
	RawSlice = Slice2(raw, 0, i, 0, 0, Nx, 1, 2, Nz, 1)
	Slice2(EnergyARCL, 0, i, 0, 0, Nx, 1, 2, Ne, 1) = Warp(RawSlice, icol, Row)
	
	// Result("\n" + i)
}

Image EnergyMomentumCL := RealImage("Energy-Momentum CL", 4, Nm, Ne, Nphi)
EnergyMomentumCL.ImageSetDimensionScale(0, Kmax_um/(Nm/2))
EnergyMomentumCL.ImageSetDimensionOrigin(0, -Kmax_um)
EnergyMomentumCL.ImageSetDimensionUnitString(0, "/um")
EnergyMomentumCL.ImageSetDimensionScale(1, -(EnergyMax_eV - EnergyMin_eV)/Ne)
EnergyMomentumCL.ImageSetDimensionOrigin(1, EnergyMax_eV)
EnergyMomentumCL.ImageSetDimensionUnitString(1, "eV")
EnergyMomentumCL.ImageSetDimensionScale(2, 180/Nphi)
EnergyMomentumCL.ImageSetDimensionOrigin(2, 0)
EnergyMomentumCL.ImageSetDimensionUnitString(2, "deg")
//EnergyMomentumCL.showImage()

Image ARCLSlice := RealImage("Current ARCL slice", 4, Nx, Ny)
//ARCLSlice.showImage()
Image EnergyMomentumSlice := RealImage("Current Energy-Momentum CL slice", 4, Nm, Nphi)
//EnergyMomentumSlice.showImage()

Image K_um := RealImage("Sample plane momentum map, K", 4, Nm, Nphi)
Image Theta := RealImage("Theta [deg]", 4, Nm, Nphi)
Image Phi := RealImage("Phi", 4, Nm, Nphi)

Image XX := RealImage("XX", 4, Nm, Nphi)
Image YY := RealImage("ZZ", 4, Nm, Nphi)
Image TT := RealImage("TT", 4, Nm, Nphi)
//XX.showImage()


//K_um = (iradius/(Nm/2))*Kmax_um
K_um = (icol/Nm - 0.5)*2*Kmax_um
//K_um.showImage()

Number Wavelength_um

number cutoff = (3.45 - energymin_ev)*ne/(energymax_ev - energymin_ev)

for (number i = 0; i < Ne; i++)
{
	// Slice2(img, Start X, Start Y, Start Z, Dim 2, steps 2, stride 2, Dim 2, Steps 2, stride 2)
	ARCLSlice = Slice2(EnergyARCL, 0, 0, i, 0, Nx, 1, 1, Ny, 1)
	
	Wavelength_um = 1.2398/(EnergyMin_eV + (EnergyMax_eV - EnergyMin_eV)*i/Ne)
	Theta = asin(K_um*Wavelength_um/(2*pi())) // Use for Momentum mapping
	Phi = (irow/Nphi)*pi()

	XX = cos(Phi)*(Theta/(AngleSampling_deg*pi()/180)) + Nx/2
	YY = sin(Phi)*(Theta/(AngleSampling_deg*pi()/180)) + Ny/2
	
	EnergyMomentumSlice = Warp(ARCLSlice, XX, YY) 
		if (i > cutoff)
		{
			EnergyMomentumSlice = Energymomentumslice*0
		}
	Slice2(EnergyMomentumCL, 0, (Ne - i - 1), 0, 0, Nm, 1, 2, Nphi, 1) = EnergyMomentumSlice
	Result("\nRemapping slice: " + i + " of " + Ne)
}


EnergyMomentumCL.showImage()