Convert real space in reciprocal space

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3     // Use modern global access method and strict wave access.

function SetScaleMap(srcWave)
    String srcWave
    variable tilt=0,phi=0
    Make/O/N=(DimSize($srcWave,1),DimSize($srcWave,2)) WXX,WYY
    SetScale/P x,DimOffset($srcWave,1),DimDelta($srcWave,1), WXX,WYY
    SetScale/P y,DimOffset($srcWave,2),DimDelta($srcWave,2), WXX,WYY
    Make/O/N=(DimSize($srcWave,0),DimSize($srcWave,1),DimSize($srcWave,2)) Test1,Test2
    SetScale/P x,DimOffset($srcWave,0),DimDelta($srcWave,0), Test1,Test2
    SetScale/I y,wavemin(wxx),wavemax(wxx), Test1,Test2
    SetScale/I z,wavemin(wyy),wavemax(wyy), Test1,Test2
    wave w=$srcWave
    multithread Test2[][][]=w[scaleToIndex(w,x,0)][scaleToIndex(w,180/pi*asin((sin(tilt*pi/180)*sqrt((0.512*sqrt(x))^2-y^2-z^2)-cos(tilt*pi/180)*(cos(phi*pi/180)*y+sin(phi*pi/180)*z))/(0.512*sqrt(x))),1)][scaletoIndex(w,180/pi*atan((sin(phi*pi/180)*y-cos(phi*pi/180)*z)/(sin(tilt*pi/180)*cos(phi*pi/180)*y+sin(tilt*pi/180)*sin(phi*pi/180)*z+cos(tilt*pi/180)*sqrt((0.512*sqrt(x))^2-y^2-z^2))),2)]


function TransitionMatrix(axe,alpha,theta,tilt,phi)
    variable axe,alpha,theta,tilt,phi
    Variable l_11=cos(phi*pi/180)*cos(tilt*pi/180)
    Variable l_12=cos(phi*pi/180)*sin(theta*pi/180)*sin(tilt*pi/180)-sin(phi*pi/180)*cos(theta*pi/180)
    Variable l_13=cos(phi*pi/180)*cos(theta*pi/180)*sin(tilt*pi/180)+sin(phi*pi/180)*sin(theta*pi/180)
    Variable l_21=sin(phi*pi/180)*cos(tilt*pi/180)
    Variable l_22=sin(phi*pi/180)*sin(theta*pi/180)*sin(tilt*pi/180)+cos(phi*pi/180)*cos(theta*pi/180)
    Variable l_23=sin(phi*pi/180)*cos(theta*pi/180)*sin(tilt*pi/180)-cos(phi*pi/180)*sin(theta*pi/180)
    Variable l_31=-sin(tilt*pi/180)
    Variable l_32=sin(theta*pi/180)*cos(tilt*pi/180)
    Variable l_33=cos(theta*pi/180)*cos(tilt*pi/180)
    variable kx=(-l_11*sin(alpha*pi/180)+l_12*0+l_13*cos(alpha*pi/180))
    variable ky=(-l_21*sin(alpha*pi/180)+l_22*0+l_23*cos(alpha*pi/180))
    variable kz=(-l_31*sin(alpha*pi/180)+l_32*0+l_33*cos(alpha*pi/180))
        return kx
        return ky



Dear IgorPro users,

I am looking for a procedure to convert real space to reciprocal space. In practice, I have data in the (xyz) layer and I will transfer each cell into a new 3D wave (kxkykz) but I have a problem. In fact, my data is ARPES data (photoemission), the new 3D wave is bigger than the original 3D wave because of the space transformation, so when I assign the new 3D wave, some data is missing in the original wave and igorpro returns me an error.

I have a solution using the interp3D function but I want to do the same operation without this function because Interp3D takes time. Do you have any ideas?


Thanks in advance


I attach my code

Can you post the analytical equations, for example in MathJax format. In the meantime, you could likely speed up the code by converting some of the repeating functions to variables at the start. For example ap = alpha*pi/180.

I have restructured your code a bit to make it more readable (for me; see below). First, the suggestion from jj is a good one. You may want to replace repeated calculations by a bit more memory usage by putting trigonometric expressions into variables as much as possible. Note that I haven't actually done that in below code, since the variables are recalculated for each call of convertY, convertZ.

Regarding your question, first you have to note that you are loosing fidelity by using the expression of test2, since you only select integer data points with this approach and not interpolated ones like with interp3D. But if this is good enough for you you can just hack your way to the result by using

#pragma rtGlobals=1

This setting will just clip any data reference to srcWave should the index be out of bounds (i.e., it suppresses the error you are getting without improving the situation). A bit more sophisticated approach would be to use a lookup function instead of direct wave referencing (e.g., something like Test2 = myLookUpFunc(w, tilt, phi, x,y,z)), which checks each index against the wave boundaries of srcWave, and may return NaN if out of bounds, and thus prevents the error from happening in the frist place. This should give you similar speed and also fills the out-of-bounds areas with a proper value.

Also, you don't seem to do much with the third dimension (x), which makes this essentially a 2D problem. I wonder if it would be faster to just use interp2D on each x plane.

function SetScaleMap(srcWave)
    String srcWave
    wave w=$srcWave
    variable tilt=0, phi=0
    variable sx = DimSize($srcWave,0), sy = DimSize($srcWave,1), sz = DimSize($srcWave,2)
    variable ox = DimOffset($srcWave,0), oy = DimOffset($srcWave,1), oz = DimOffset($srcWave,2)
    variable dx = DimDelta($srcWave,0), dy = DimDelta($srcWave,1), dz = DimDelta($srcWave,2)   
    Make/O/N=(sy,sz) WXX,WYY
    SetScale/P x,oy,dy, WXX,WYY
    SetScale/P y,oz,dz, WXX,WYY
    wxx = 0.512*sqrt((sx-1)*dx+ox) * TransitionMatrix(0,x,y,tilt,phi)
    wyy = 0.512*sqrt((sx-1)*dx+ox) * TransitionMatrix(1,x,y,tilt,phi)
    Make/O/N=(sx,sy,sz) Test1,Test2
    SetScale/P x,ox,dx, Test1,Test2
    SetScale/I y,wavemin(wxx),wavemax(wxx), Test1,Test2
    SetScale/I z,wavemin(wyy),wavemax(wyy), Test1,Test2
    Test1 = interp3D(w,x,convertY(tilt,phi,x,y,z),convertZ(tilt,phi,x,y,z))
    multithread Test2[][][] = w[scaleToIndex(w,x,0)][scaleToIndex(w,convertY(tilt,phi,x,y,z),1)][scaletoIndex(w,convertZ(tilt,phi,x,y,z),2)]


threadsafe function convertY(tilt,phi,x,y,z)
    variable tilt,phi,x,y,z

    variable ct = cos(tilt*pi/180), cp = cos(phi*pi/180)
    variable st = sin(tilt*pi/180), sp = sin(phi*pi/180)

    return 180/pi*asin((st*sqrt((0.512*sqrt(x))^2-y^2-z^2)-ct*(cp*y+sp*z))/(0.512*sqrt(x)))

threadsafe function convertZ(tilt,phi,x,y,z)
    variable tilt,phi,x,y,z

    variable ct = cos(tilt*pi/180), cp = cos(phi*pi/180)
    variable st = sin(tilt*pi/180), sp = sin(phi*pi/180)

    return 180/pi*atan((sp*y-cp*z)/(st*cp*y+st*sp*z+ct*sqrt((0.512*sqrt(x))^2-y^2-z^2)))

function TransitionMatrix(axe,alpha,theta,tilt,phi)
    variable axe,alpha,theta,tilt,phi
    Variable l_11 = cos(phi*pi/180)*cos(tilt*pi/180)
    Variable l_12 = cos(phi*pi/180)*sin(theta*pi/180)*sin(tilt*pi/180) - sin(phi*pi/180)*cos(theta*pi/180)
    Variable l_13 = cos(phi*pi/180)*cos(theta*pi/180)*sin(tilt*pi/180) + sin(phi*pi/180)*sin(theta*pi/180)
    Variable l_21 = sin(phi*pi/180)*cos(tilt*pi/180)
    Variable l_22 = sin(phi*pi/180)*sin(theta*pi/180)*sin(tilt*pi/180) + cos(phi*pi/180)*cos(theta*pi/180)
    Variable l_23 = sin(phi*pi/180)*cos(theta*pi/180)*sin(tilt*pi/180) - cos(phi*pi/180)*sin(theta*pi/180)
    Variable l_31 = -sin(tilt*pi/180)
    Variable l_32 = sin(theta*pi/180)*cos(tilt*pi/180)
    Variable l_33 = cos(theta*pi/180)*cos(tilt*pi/180)
    variable kx = (-l_11*sin(alpha*pi/180)+l_12*0+l_13*cos(alpha*pi/180))
    variable ky = (-l_21*sin(alpha*pi/180)+l_22*0+l_23*cos(alpha*pi/180))
    variable kz = (-l_31*sin(alpha*pi/180)+l_32*0+l_33*cos(alpha*pi/180))
        return kx
        return ky


I cannot work out the math from this. Equations would be helpful. If this could be somehow calculated using MatrixOP, it usually is much faster. When I needed to convert between real and reciprocal space, FFT and IFFT have been helpful - and much faster. 

Some efficiency improvements which will not solve the issue, but may speed it up enough:

You could use all angles in radians directly and avoid all conversions from degrees. Would make everything faster and more readable. 

In TransitionMatrix function you can precalculate all those cos(phi*pi/180) and others.

kz calculated in TransitionMatrix is never even used, so there is bunch of lines which can be removed to avoid spending time on them. 

Interp3D is multithreadable (at least in IP9), can you simply add multithread in front of the Test1= line?  (would need to declare some functions threadsafe, but I think they are). 

Thanks for the answers.

Actually, my real space is not exactly a real space, it is an angular and energetic space. Due to conservation laws, this space can be converted to a reciprocal (kx,ky,kz) space ( FFT and IFFT are therefore not appropriate here, sorry for the misuse of language.

During a measurement, we have to change the orientation of the sample, which creates a distortion of the data. It is possible to correct this distortion with a change basis ( So I want to implement this idea in an igorpro procedure. 

The first solution uses an interpolation of the data along a certain path (this is a 2D problem, not 3D, you are right). The second solution is to define a new 3D wave (kx,ky,kz) and to assign each cell of the original wave to a cell of the new wave via a search function. It is this second way that causes me problems. I hope to have clarified a little my project. I enclose a capture of the data before and after correction made with the interpolation function.




As mentioned, you have at least four options to speed up things:

- use interp3D and try to optimize your code for speed

- try to use interp2D and see how far this gets you

- don't use interpolation and just ignore the error

- don't use interpolation and write an assignment function which handles the data range properly

Whichever you choose, people in the forum can help you to get the code working. Of course, there is the fifth option to leave things as is and just let the code running (you didn't mention how long you have to wait).

Yes, I tried to use interp2D in the for loop but it takes longer than interp3D. The interp3D takes approximately twenty seconds to process the data layers (600x600x600) (with the multithreaded option on IP8). This is not so long and it should be able to improve with the improvement of my code.
Now I want to work on the fourth possibility but I have no idea how to start the lookup function.

Here this should work (not tested):

Threadsafe Function indexLookUp(w,tilt,phi,wx,wy,wz)
    Wave w
    variable tilt,phi,wx,wy,wz
    variable wp = scaleToIndex(w,wx,0)
    variable wq = scaleToIndex(w,convertY(tilt,phi,wx,wy,wz),1)
    variable wr = scaletoIndex(w,convertZ(tilt,phi,wx,wy,wz),2)
    if (wq < 0 || wq > DimSize(w,1)-1 || wr < 0 || wr > DimSize(w,2)-1)
        return NaN
        return w[wp][wq][wr]

then do:

multithread Test2 = indexLookUp(w,tilt,phi,wx,wy,wz)

This basically prevents access to srcwave if the index is out of bounds.

Yes, it works.
Now I want to do the search function. The idea is to create a wave wxx,wyy from the wave scale data and create a function to search for that point on new3Dwave. If true, the new cell is assigned from the data, otherwise Nan is assigned to the cell. Currently I have four forloops to browse my wave, but it's not very efficient.

It is a bit difficult to imagine what you want to do from your description. Could you post some code?