ARPES: E(deg) to E(K)
Tue, 12/29/2020 - 02:05 pm
I usually use the code that I'll share below to obtain the final ARPES image.
The problem is when the data contains many points it takes a lot of time to do it.
I'll be grateful if you could help with ideas.
Here is the code:
#pragma rtGlobals=1 // Use modern global access method.
#include <XYZtoMatrix>
Function/S ARPES(DataWave)
Wave DataWave;
// Duplicate the wave for backup, and rename it to RealSpace
Duplicate DataWave RealSpace;
// Determine the size of two dimensional wave
Variable NumberOfRow;
Variable NumberOfColumn;
NumberOfRow = DimSize(RealSpace, 0);
NumberOfColumn = DimSize(RealSpace, 1);
// Rescale the 'Scaled dimension index' of YScale into rad
Variable YScaleOffset;
YScaleOffset = DimOffset(RealSpace, 1);
YScaleOffset = Pi * YScaleOffset / 180;
Variable YScaleDelta;
YScaleDelta = DimDelta(RealSpace, 1);
YScaleDelta = Pi * YScaleDelta / 180;
SetScale/P y YScaleOffset, YScaleDelta, "rad", RealSpace;
// Extract the rad values to a new wave
Make/N=(NumberOfColumn)/D RadValue;
RadValue = YScaleOffset + p*YScaleDelta;
// Determine the DimOffset and DimDelta of Column
Variable XScaleOffset;
Variable XScaleDelta;
XScaleOffset = DimOffset(RealSpace, 0);
XScaleDelta = Dimdelta(RealSpace, 0);
// Length of k_space vectors
Variable LengthOfKSpaceWave;
LengthOfKSpaceWave = NumberOfRow * NumberOfColumn;
// Create three waves of length LengthOfKSpaceWave for E_B, K_parallel, and Counts
Make/N=(LengthOfKSpaceWave)/D E_B;
Make/N=(LengthOfKSpaceWave)/D K_parallel;
Make/N=(LengthOfKSpaceWave)/D Counts;
Variable ii;
Variable jj;
Variable a;
Variable E_F;
E_F = 0;
Prompt E_F, " Enter Fermi enegy"
//Prompt Field_Oe, "Applied magnetic field (Oe):"
DoPrompt "Please enter the Fermi enegy (eV)", E_F;
For(jj=0; jj<(NumberOfRow); jj+=1)
For(ii=0; ii<(NumberOfColumn); ii+=1)
a = (jj*NumberOfColumn) + ii;
K_parallel[a] = 0.512*Sqrt(XScaleOffset + XScaleDelta*jj)*Sin(RadValue[ii]);
E_B[a] = (XScaleOffset + XScaleDelta*jj) - E_F;
Counts[a] = RealSpace[jj][ii];
EndFor
//String NewCountWaveName = "Counts_"+ num2str(ii);
//String NewK_parallelWaveName = "K_parallel_" + num2str(ii);
//Rename K_parallel, $NewK_parallelWaveName;
//Rename Counts, $NewCountWaveName;
EndFor
//KillWaves NewK_parallelWaveName, NewCountWaveName;
KillWaves RadValue;
KillWaves RealSpace;
Display K_parallel vs E_B;
ModifyGraph mode = 3, Marker = 19;
ModifyGraph zColor(K_parallel)={Counts,*,*,Grays,0};
Label left "K_parallel";DelayUpdate
Label bottom "E_B (eV)"
End
Can you edit this code and put it inside the Code sniped block? It will be much easier to read.
My best guess is that most time is spent inside the two embedded for loops. If this can be written as Igor wave line, it could be called multithreaded. If it cannot be written as one line, then may be as threadsafe function and called multithreaded. You would get much higher performance in either case. I do not think rest of the code takes much time.
There are few tricks you can surely make easier - there is no need for K_parallel, E_B, and Counts to be 1D vectors in that calculation, you can have them as 2D matrix, assign values in each [p][q] and convert to vector outside the loop by smart use of redimension. That will make it possible to write each as simple wave assigmnent (single line).
I think these: (XScaleOffset + XScaleDelta*jj) are simply X value which Igor gives you already? If not, why not make it somehow...
Use free waves, they seem to be faster to deal with and no need to kill them eventually.
December 29, 2020 at 02:22 pm - Permalink
#include <XYZtoMatrix>
Function/S ARPES(DataWave)
Wave DataWave;
// Duplicate the wave for backup, and rename it to RealSpace
Duplicate DataWave RealSpace;
// Determine the size of two dimensional wave
Variable NumberOfRow;
Variable NumberOfColumn;
NumberOfRow = DimSize(RealSpace, 0);
NumberOfColumn = DimSize(RealSpace, 1);
// Rescale the 'Scaled dimension index' of YScale into rad
Variable YScaleOffset;
YScaleOffset = DimOffset(RealSpace, 1);
YScaleOffset = Pi * YScaleOffset / 180;
Variable YScaleDelta;
YScaleDelta = DimDelta(RealSpace, 1);
YScaleDelta = Pi * YScaleDelta / 180;
SetScale/P y YScaleOffset, YScaleDelta, "rad", RealSpace;
// Extract the rad values to a new wave
Make/N=(NumberOfColumn)/D RadValue;
RadValue = YScaleOffset + p*YScaleDelta;
// Determine the DimOffset and DimDelta of Column
Variable XScaleOffset;
Variable XScaleDelta;
XScaleOffset = DimOffset(RealSpace, 0);
XScaleDelta = Dimdelta(RealSpace, 0);
// Length of k_space vectors
Variable LengthOfKSpaceWave;
LengthOfKSpaceWave = NumberOfRow * NumberOfColumn;
// Create three waves of length LengthOfKSpaceWave for E_B, K_parallel, and Counts
Make/N=(LengthOfKSpaceWave)/D E_B;
Make/N=(LengthOfKSpaceWave)/D K_parallel;
Make/N=(LengthOfKSpaceWave)/D Counts;
Variable ii;
Variable jj;
Variable a;
Variable E_F;
E_F = 0;
Prompt E_F, " Enter Fermi enegy"
//Prompt Field_Oe, "Applied magnetic field (Oe):"
DoPrompt "Please enter the Fermi enegy (eV)", E_F;
For(jj=0; jj<(NumberOfRow); jj+=1)
For(ii=0; ii<(NumberOfColumn); ii+=1)
a = (jj*NumberOfColumn) + ii;
K_parallel[a] = 0.512*Sqrt(XScaleOffset + XScaleDelta*jj)*Sin(RadValue[ii]);
E_B[a] = (XScaleOffset + XScaleDelta*jj) - E_F;
Counts[a] = RealSpace[jj][ii];
EndFor
//String NewCountWaveName = "Counts_"+ num2str(ii);
//String NewK_parallelWaveName = "K_parallel_" + num2str(ii);
//Rename K_parallel, $NewK_parallelWaveName;
//Rename Counts, $NewCountWaveName;
EndFor
//KillWaves NewK_parallelWaveName, NewCountWaveName;
KillWaves RadValue;
KillWaves RealSpace;
Display K_parallel vs E_B;
ModifyGraph mode = 3, Marker = 19;
ModifyGraph zColor(K_parallel)={Counts,*,*,Grays,0};
Label left "K_parallel";DelayUpdate
Label bottom "E_B (eV)"
End
The first part was easy.
Andy
December 29, 2020 at 02:29 pm - Permalink
And you don't need the ending semicolons.
December 29, 2020 at 03:23 pm - Permalink
In reply to Can you edit this code and… by ilavsky
#include <XYZtoMatrix>
Function/S ARPES(DataWave)
Wave DataWave;
// Duplicate the wave for backup, and rename it to RealSpace
Duplicate DataWave RealSpace;
// Determine the size of two dimensional wave
Variable NumberOfRow;
Variable NumberOfColumn;
NumberOfRow = DimSize(RealSpace, 0);
NumberOfColumn = DimSize(RealSpace, 1);
// Rescale the 'Scaled dimension index' of YScale into rad
Variable YScaleOffset;
YScaleOffset = DimOffset(RealSpace, 1);
YScaleOffset = Pi * YScaleOffset / 180;
Variable YScaleDelta;
YScaleDelta = DimDelta(RealSpace, 1);
YScaleDelta = Pi * YScaleDelta / 180;
SetScale/P y YScaleOffset, YScaleDelta, "rad", RealSpace;
// Extract the rad values to a new wave
Make/N=(NumberOfColumn)/D RadValue;
RadValue = YScaleOffset + p*YScaleDelta;
// Determine the DimOffset and DimDelta of Column
Variable XScaleOffset;
Variable XScaleDelta;
XScaleOffset = DimOffset(RealSpace, 0);
XScaleDelta = Dimdelta(RealSpace, 0);
// Length of k_space vectors
Variable LengthOfKSpaceWave;
LengthOfKSpaceWave = NumberOfRow * NumberOfColumn;
// Create three waves of length LengthOfKSpaceWave for E_B, K_parallel, and Counts
Make/N=(LengthOfKSpaceWave)/D E_B;
Make/N=(LengthOfKSpaceWave)/D K_parallel;
Make/N=(LengthOfKSpaceWave)/D Counts;
Variable ii;
Variable jj;
Variable a;
Variable E_F;
E_F = 0;
Prompt E_F, " Enter Fermi enegy"
//Prompt Field_Oe, "Applied magnetic field (Oe):"
DoPrompt "Please enter the Fermi enegy (eV)", E_F;
For(jj=0; jj<(NumberOfRow); jj+=1)
For(ii=0; ii<(NumberOfColumn); ii+=1)
a = (jj*NumberOfColumn) + ii;
K_parallel[a] = 0.512*Sqrt(XScaleOffset + XScaleDelta*jj)*Sin(RadValue[ii]);
E_B[a] = (XScaleOffset + XScaleDelta*jj) - E_F;
Counts[a] = RealSpace[jj][ii];
EndFor
//String NewCountWaveName = "Counts_"+ num2str(ii);
//String NewK_parallelWaveName = "K_parallel_" + num2str(ii);
//Rename K_parallel, $NewK_parallelWaveName;
//Rename Counts, $NewCountWaveName;
EndFor
//KillWaves NewK_parallelWaveName, NewCountWaveName;
KillWaves RadValue;
KillWaves RealSpace;
Display K_parallel vs E_B;
ModifyGraph mode = 3, Marker = 19;
ModifyGraph zColor(K_parallel)={Counts,*,*,Grays,0};
Label left "K_parallel";DelayUpdate
Label bottom "E_B (eV)"
End
December 29, 2020 at 03:43 pm - Permalink
In reply to Can you edit this code and… by ilavsky
Exactly, that's the point! how to make it 2D matrix directly
I didn't understand how to convert into vector
December 29, 2020 at 04:21 pm - Permalink
This is my wild guess without testing... This is for K_parallel only, the E_B is similar. But I think this like this:
Duplicate/O DataWave, K_parallel - you will have 2D K_parallel with the same x and y scaling as the data.
Then I think your two nested loops do simply
K_parallel = 0.512*sqrt(x)*sin(y)
I may be wrong, but I think this is what you are doing. And if this takes too long, do this in with multithread keyword in front of the equation and it will run on all cores you have:
multithread K_parallel = 0.512*sqrt(x)*sin(y)
After that you do
Redimension/N=1 K_parallel
and it will convert 2D K_parallel into 1D K_parallel (which is anyway how Igor orders the data internally, check the manual for how this is done, it is well explained there in About Waves I think). If the order of Columns/rows in 1D output is wrong, you may need to swap the K_parallel columns/rows first using MatrixOp operation.
December 29, 2020 at 04:52 pm - Permalink
You can use interp2D effectively to convert angles to k values of an ARPES map (2D->2D). Here is the rusty code from my university days which does that. It may not work out of the box (I haven't looked into your code what your data looks like), but should give you an idea. It does not get rid of the nested loop, but I guess you can optimize this code for speed in many places.
Wave inwave
Variable rows,columns,xdelta,xoffset,ydelta,yoffset
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,prec = 2
Emax = xoffset + xdelta*(rows-1)
kmin = 0.512*sqrt(Emax)*sin(pi/180*(yoffset))
kmax = 0.512*sqrt(Emax)*sin(pi/180*(yoffset+(columns-1)*ydelta))
String NewName = NameofWave(inwave)+"_k"
Make/O/N=(rows,prec*columns) $NewName
Wave outwave = $NewName
SetScale/P x xoffset,xdelta, WaveUnits(inwave,0), outwave
SetScale/I y kmin,kmax, "k [Å-1]", outwave
kdelta = DimDelta(outwave,1)
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Variable i,j,IPx,IPy, value
for (i = 0; i < rows; i += 1)
for (j = 0; j < prec*columns; j += 1)
IPx = xoffset + i*xdelta
IPy = 180/pi*asin( (kmin+j*kdelta) / ( 0.512*sqrt(IPx) ) )
value = interp2D(inwave, IPx, IPy)
if (value < 0 || value > 0)
outwave[i][j] = value
endif
endfor
endfor
End
December 29, 2020 at 05:34 pm - Permalink