ARPES: E(deg) to E(k) with bias voltage

I have recently seen a code written to convert the E(deg) to E(k). IgorExchange
Now I would like to change the code to accommodate the biasing effect. The new equation for E(deg) to E(k||) with biasing voltage will look like in the attached file. Although theta is function of E_kin but it is similar to common K|| conversion equation.
I have tried to write as follows: mainly my aim is to do all this things in single macro. Please let me know does it possible.

Here is the code:
Function AngleToK_withbias(inwave) /// needs angle (y-axis) vs energy (x-axis)
Wave inwave

String newname = NameofWave(inwave)+"_kwb"
Duplicate/O inwave, $newname
Wave outwave = $newname

Variable rows,columns,xdelta,xoffset,ydelta,yoffset // inwave parameters
rows = DimSize(inwave,0)
columns = DimSize(inwave,1)
xdelta = DimDelta(inwave,0)
xoffset = DimOffset(inwave,0)
ydelta = DimDelta(inwave,1)
yoffset = DimOffset(inwave,1)

Variable kmin,kmax,kdelta,Emax,Emin
Emax = xoffset + xdelta*(rows-1)
Variable theta_f,theta_i
Variable U=5
Variable e= 1.602176565*10^-19
theta_i=pi/180*(yoffset)
Emin=xoffset
theta_f=theta_i+atan((e*U*(1-sin(theta_i)/theta_i))/(Emin-(e*U*(1-sin(theta_i)/theta_i))))^0.5
kmin = 0.512*sqrt(Emax)*sin(theta_f) // calculate the k boundaries (i.e., k at highest Ekin)
theta_i=pi/180*(yoffset+(columns-1)*ydelta)
theta_f=theta_i+atan((e*U*(1-sin(theta_i)/theta_i))/(Emin-(e*U*(1-sin(theta_i)/theta_i))))^0.5
kmax = 0.512*sqrt(Emax)*sin(theta_f)
SetScale/I y kmin,kmax,"Ang^-1", outwave // scale the y axis

outwave = interp2D(inwave, x, (180/pi)*asin(y/ (0.512*sqrt(x)))) // recalculate to k
outwave = (NumType(outwave)==2) ? 0 : outwave // replace NaNs (optional)
End

If you have any idea pls let me know.
formula.pdf
I have written a function to apply the angle correction only. It is rather difficult to integrate also the k conversion in one step since the bias correction is rather convoluted. But it is of course very easy to apply the k vector calculation directly afterwards in the same function. I use the same conventions as for my k conversion code, i.e., x is energy and y is angle, just to be consistent. By the way I didn't use your posted formula but equation A6 from the appendix of the paper (Phys. Rev. B 77, 085425 (2008)), since A7 seems to be a simplified case. Anyway, if you like to use equation A7 just replace the simple code of the AngleFormula function. Here's the code (I also attach the procedure file):
Function BiasAngle(inwave,U,tilt)
    Wave inwave
    Variable U, tilt // U in eV, tilt in deg
   
    String newname = NameofWave(inwave)+"_b"
    Duplicate/O inwave, $newname
    Wave outwave = $newname

    Variable rows,columns,xdelta,xoffset,ydelta,yoffset     // inwave parameters
    rows    = DimSize(inwave,0)
    columns = DimSize(inwave,1)
    xoffset = DimOffset(inwave,0)
    ydelta  = DimDelta(inwave,1)
    yoffset = DimOffset(inwave,1)
   
    SetScale/I  y yoffset+AngleFormula(U,xoffset,tilt), yoffset+ydelta*(columns-1)+AngleFormula(U,xoffset,tilt), outwave    // scale the y axis to the extreme vaules
   
    outwave = interp2D(inwave, x, y - AngleFormula(U,x,tilt))   // recalculate to corrected angles
    outwave = (NumType(outwave)==2) ? 0 : outwave           // replace NaNs (optional)
End

Function AngleFormula(U,Ekin,tilt)  // output: correction in deg
    Variable U, Ekin, tilt  // U in eV, Ekin in eV, tilt in deg
    Variable radT = tilt*Pi/180
    Variable inbrackets=1-sqrt(1-radT^2)/2-asin(radT)/(2*radT)
    return asin(sqrt(U/Ekin*inbrackets))/(Pi/180)
End
BiasAngle.ipf
You're welcome. These kind of problems help to learn new stuff in Igor, so if it's not too difficult (i.e., require several hours of work) I can only benefit. ;)
Actually I discovered a little mistake in the code. The lower limit of the scaling was wrong so the data gets cut off. Also, I though it would be useful to include the case of when the tilt angle gets negative (i.e., when you have tilted the sample in the other direction but don't want to invert the data). U should be positive (although you can insert a negative value since I use the absolute).

Anyway, here's the new code:
Function BiasAngle(inwave,U,tilt)
    Wave inwave
    Variable U, tilt // U in eV, tilt in deg
   
    String newname = NameofWave(inwave)+"_b"
    Duplicate/O inwave, $newname
    Wave outwave = $newname

    Variable rows,columns,xdelta,xoffset,ydelta,yoffset     // inwave parameters
    rows    = DimSize(inwave,0)
    columns = DimSize(inwave,1)
    xdelta  = DimDelta(inwave,0)
    xoffset = DimOffset(inwave,0)
    ydelta  = DimDelta(inwave,1)
    yoffset = DimOffset(inwave,1)
   
    Variable MinE = xoffset + (tilt<0)*xdelta*(rows-1), MaxE = xoffset + (tilt>=0)*xdelta*(rows-1)
    SetScale/I y yoffset+sign(tilt)*AngleFormula(U,MaxE,tilt), yoffset+ydelta*(columns-1)+sign(tilt)*AngleFormula(U,MinE,tilt), outwave // scale the y axis to the extreme vaules
   
    outwave = interp2D(inwave, x, y - sign(tilt)*AngleFormula(U,x,tilt))    // recalculate to corrected angles
    outwave = (NumType(outwave)==2) ? 0 : outwave           // replace NaNs (optional)
End

Function AngleFormula(U,Ekin,tilt)  // output: correction in deg
    Variable U, Ekin, tilt  // U in eV, Ekin in eV, tilt in deg
    Variable radT = tilt*Pi/180
    Variable inbrackets=1-sqrt(1-radT^2)/2-asin(radT)/(2*radT)
    return asin(sqrt(abs(U)/Ekin*inbrackets))/(Pi/180)
End
As a general discussion, I would like to share with you the following: Actually, data has been taken in different emission angle to probe the more k-space. Later, they are stitched togather to get single image which will have say -60 to +60 deg range. So, I thing one may not need to feed the tilt value as input (in my case, so I will modify the code). But it is really great to have versatile program! Thanks for modified code that considers the negative angle as well.
Of course you need the tilt angle, since otherwise the whole calculation is meaningless. I guess you apply the angle shift to the tilted spectra anyway to stitch the data together. So just use the procedure (which automatically includes the angle shift) before the stitching step to correct the data. I can't imagine how you would do it after everything has been stitched, as you can't tell apart the parts of the spectrum where the sample was tilted and where not. Also I recommend to measure some overlap since the data gets skewed a bit by the bias correction.
If one measure ARPES data at normal emission mode with bias voltage, then bias effect is present but there is no tilt!
So, in that case what ever angle we have in the image in y-axis, that will be the value one should feed in to the equation as theta_m. But, when data is taken with sample is rotated then tilt is must to have correction. This is also necessary just because the y-scale is not corrected. So, when one stitch the image taken at different tilt angle and get single image with proper scaling, then angle value (y-scale) is just enough and that is the value should be feed to the equation as theta_m.
I see. I was under the impression that tilting the surface with respect to the analyzer plane leads to a different geometry of you E-field. But if that plays no role, then I guess in can be done afterwards.