4Misc_StartP#4Platform@9VersionCheck xHH@Rg(HHdh xHH@Rg(HHdh x HH@Rg(HHdh ^Graph*xpWDashSettings#  !4 4 4 4 4 4 whomeCCd Macintosh HD:Users:lingyuan:Desktop:cRecentWindows/Auxiliaryd3dErrors.ihfHaldaneHaldane Model 4Misc_EndP#4XOPState_StartP#4XOPState_EndP#\lA..C'a绦W%;!K*d TXET0RGIhHaldaneHaldane#pragma TextEncoding = "UTF-8" #pragma rtGlobals=3 // Use modern global access method and strict wave access #pragma DefaultTab={3,20,4} // Set default tab width in Igor Pro 9 and later /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// //****************************** Haldane Model ******************************************** /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// //Ref: Duncan Haldane Phys. Rev. Lett. 61 (1988) 2015-2018 //Ref: https://ncatlab.org/nlab/show/Haldane+model#DTC //Ref: https://topocondmat.org/w4_haldane/haldane_model.html /////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////// //************************* Build Haldane Hamiltonian ************************************* /////////////////////////////////////////////////////////////////////////////////////////// Function/Wave HaldaneModel(t1,t2,M,phi,kx,ky) variable t1 variable t2 variable M variable phi variable kx variable ky //** Lattice vector of Haldane Hamiltonian variable a = 1 wn("a1",{0,a}) wn("a2",{-sqrt(3)*a/2,-a/2}) wn("a3",{sqrt(3)*a/2,-a/2}) wave a1 = $"a1" wave a2 = $"a2" wave a3 = $"a3" duplicate/o a1 b1 duplicate/o a1 b2 duplicate/o a1 b3 b1 = a2-a3 b2 = a3-a1 b3 = a1-a2 variable e,dx,dy,dz //** Pointer String for called Functions string ps = "ps_haldane" string result_C = "result_C" //** Create Haldane Hamiltonian e = 2*t2*cos(phi)*(cos(kx*b1[0]+ky*b1[1]) + cos(kx*b2[0]+ky*b2[1]) + cos(kx*b3[0]+ky*b3[1])) dx = t1*(Cos(kx*a1[0]+ky*a1[1])+Cos(kx*a2[0]+ky*a2[1])+Cos(kx*a3[0]+ky*a3[1])) dy = t1* (Sin(kx*a1[0]+ky*a1[1])+Sin(kx*a2[0]+ky*a2[1])+Sin(kx*a3[0]+ky*a3[1])) dz = M-2*t2*Sin(phi)*(Sin(kx*b1[0]+ky*b1[1])+Sin(kx*b2[0]+ky*b2[1])+Sin(kx*b3[0]+ky*b3[1])) //## Haldane Hamiltonian: H=e*sigma0+dx*sigma1+dy*sigma2+dz*sigma3; wC("e_haldane",{{cmplx(e,0)}}) wC("dx_haldane",{{cmplx(dx,0)}}) wC("dy_haldane",{{cmplx(dy,0)}}) wC("dz_haldane",{{cmplx(dz,0)}}) Wn("ps_haldane",{{0},{1},{2},{3}}) automatrixC("e_haldane;dx_haldane;dy_haldane;dz_haldane",$ps) wave/C result_Cw = $result_C duplicate/C/o result_Cw HaldaneH killwaves $result_C a1 a2 a3 b1 b2 b3 return HaldaneH end /////////////////////////////////////////////////////////////////////////////////////////// //*************************** Build Haldane Spectral Function ***************************** /////////////////////////////////////////////////////////////////////////////////////////// Function ButtonProc_HaldaneA(ctrlName) : ButtonControl String ctrlName execute "HaldaneSpectralFunctionc()" end Proc HaldaneSpectralFunctionc(t1,t2,M,phi,selfimg,numxy) variable t1 = 1 variable t2 = 0.2 variable M = 0.2 variable phi = 0.2*pi variable selfimg = 0.05 variable numxy = 100 prompt selfimg,"Im(Σ)" prompt numxy,"N x N in (kx,ky) space" HaldaneSpectralFunction(t1,t2,M,phi,selfimg,numxy) end Function HaldaneSpectralFunction(t1,t2,M,phi,selfimg,numxy) variable t1 //= 1 variable t2 //= 0 variable M //= 0.2 variable phi //= 0.2*pi variable selfimg variable numxy //= 100 variable/G zn_cons = 2 //** Make void matrix for spectral function make/o/n=(numxy,numxy,500) modelresult_haldane setscale/i x,-3,3,"",modelresult_haldane setscale/i y,-3,3,"",modelresult_haldane setscale/i z,-4.5,4.5,"",modelresult_haldane modelresult_haldane = 0 variable roff = dimoffset(modelresult_haldane,2) variable rdelta = dimdelta(modelresult_haldane,2) //** Variable for Spectral Function variable i,j,selfreal,t selfreal = 0 variable kx variable ky //** Pointer String for called Functions string W_eigenvalues="W_eigenvalues" string HaldaneH = "HaldaneH" //** Start loop for (kx, ky) i=0 do j=0 do //** Define kx, ky by void matrix shape kx = dimoffset(modelresult_haldane,0)+i*dimdelta(modelresult_haldane,0) ky = dimoffset(modelresult_haldane,1)+j*dimdelta(modelresult_haldane,1) HaldaneModel(t1,t2,M,phi,kx,ky) wave/C HaldaneHw = $HaldaneH //** Solve eigenvalue of Haldane Hamiltonian matrixEigenV HaldaneHw wave/C n=$W_eigenvalues make/n=(dimsize(HaldaneHw,0))/o sorteigen sorteigen[0]= real(N[0]) sorteigen[1]= real(N[1]) //Make 3D Spectral Function t=0 do modelresult_haldane[i][j][]+=(-1/pi)*selfimg/(selfimg^2+((roff+r*rdelta)-selfreal-sorteigen[t])^2) t+=1 while (t < dimsize(sorteigen,0)) j+=1 while(j < dimsize(modelresult_haldane,1)) i+=1 while(i < dimsize(modelresult_haldane,0) ) //** Call smart 3D displayer duplicate/o modelresult_haldane N1_002 killwaves modelresult_haldane sorteigen HaldaneHw n execute "d3dsimu(\"N1_002\",2)" //** Modify the 3D player for this special simulation use (band structure) variable/G divcolor_cons divcolor_cons = 1 wavestats/Q $tpw() variable rangls = max(abs(V_max),abs(V_min)) ModifyImage/W=$grabwinchild(tpw()) $tpw() ctab= {-rangls,rangls,:Packages:NewColortable:dvg_seismic,0}; ColorScale/W=$grabwinchild(tpw())/C/N=textcc/X=-21.00/Y=5.00/F=0 heightPct=30,nticks=5,minor=1,tickLen=1.00,image=$tpw()//;ColorScale/W=$grabwinchild(Fexyd)/C/N=text0/F=0/A=MC/Y=-120 image=$Fexyd variable/G normornot_3dplot =2 Button launchLinecut title="Linecut",size={65,14},fSize=10,proc=ButtonProc_Cons3dplotZ end /////////////////////////////////////////////////////////////////////////////////////////// //************************** Calculate Berry Curvature ************************************ /////////////////////////////////////////////////////////////////////////////////////////// Function ButtonProc_HaldaneBC(ctrlName) : ButtonControl String ctrlName execute "HaldaneBerryCurvaturec()" end Proc HaldaneBerryCurvaturec(t1,t2,M,phi,numxy,sel) variable t1 = 1 variable t2 = 0.2 variable M = 0.2 variable phi = 0.2*pi variable numxy = 100 variable sel = 1 prompt numxy,"N x N in (kx,ky) space" prompt sel,"Select algorithm",popup,"Direct;Fukui" if (sel == 1) HaldaneBerryCurvature(t1,t2,M,phi,numxy) endif if (sel == 2) HaldaneBCFukui(t1,t2,M,phi,numxy) endif end Function HaldaneBerryCurvature(t1,t2,M,phi,numxy) variable t1 //= 1 variable t2 //= 0.2 variable M //= 0.2 variable phi //= 0.2*pi variable numxy //= 100 variable kx variable ky variable i,j Variable timerRefNum timerRefNum = StartMSTimer //** Make void matrix for spectral function make/o/n=(numxy,numxy,1) voidberry_haldane setscale/i x,-3,3,"",voidberry_haldane setscale/i y,-3,3,"",voidberry_haldane //setscale/i z,-2.5,4.5,"",voidberry_haldane voidberry_haldane = 0 //** Make void matrix for Berry curvature band #1 (lower) and #2 (higher) make/o/n=(numxy,numxy) BC_haldane1 setscale/i x,-3,3,"",BC_haldane1 setscale/i y,-3,3,"",BC_haldane1 BC_haldane1 = 0//cmplx(0,0) make/o/n=(numxy,numxy) BC_haldane2 setscale/i x,-3,3,"",BC_haldane2 setscale/i y,-3,3,"",BC_haldane2 BC_haldane2 = 0//cmplx(0,0) //** Pointer String for called Functions string W_eigenvalues="W_eigenvalues" string HaldaneH = "HaldaneH" string M_R_eigenVectors = "M_R_eigenVectors" string M_product = "M_product" //** variable for Berry curvature variable E1k, E2K //1 for low band, 2 for high band variable delta = 10^-6 variable/C Berryterm1A = cmplx(0,0) variable/C Berryterm1B = cmplx(0,0) variable/C Berryterm2A = cmplx(0,0) variable/C Berryterm2B = cmplx(0,0) //** Start loop for (kx, ky) i=0 do j=0 do //** Define kx, ky by void matrix shape kx = dimoffset(voidberry_haldane,0)+i*dimdelta(voidberry_haldane,0) ky = dimoffset(voidberry_haldane,1)+j*dimdelta(voidberry_haldane,1) HaldaneModel(t1,t2,M,phi,kx,ky) wave/C HaldaneHw = $HaldaneH //** Solve eigenvalue of Haldane Hamiltonian matrixEigenV/R HaldaneHw wave/C n=$W_eigenvalues wave/C eignevector=$M_R_eigenVectors make/n=(dimsize(HaldaneHw,0))/o sorteigen sorteigen[0]= real(N[0]) sorteigen[1]= real(N[1]) //**Sort eigenvalue and eigenvector make/o/n=(dimsize(HaldaneHw,0)) sortindex sortindex=p sort sorteigen sortindex //the value indicates the column number in M_R_eigenVectors belong to each eigenvalue labeled in sorteigen sort sorteigen sorteigen //** Calculat Berry Curvature //## E1,k = sorteigen[0] //## |u1,k> = eignevector[][sortindex[0]] //## E2,k = sorteigen[1] //## |u2,k> = eignevector[][sortindex[1]] //** get eigenvalue E1k = sorteigen[0] E2k = sorteigen[1] //** make |u1,k> //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u2k_ket u2k_ket[] = eignevector[p][sortindex[1]] // = eignevector[][sortindex[0]] //## E2,k = sorteigen[1] //## |u2,k> = eignevector[][sortindex[1]] //** get eigenvalue E1k = sorteigen[0] E2k = sorteigen[1] //** make |u1,k> //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u2k_ket u2k_ket[] = eignevector[p][sortindex[1]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u2k_ket u2k_ket[] = eignevector[p][sortindex[1]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u1kdx_ket u1kdx_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u2kdx_ket u2kdx_ket[] = eignevector[p][sortindex[1]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u1kdy_ket u1kdy_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u2kdy_ket u2kdy_ket[] = eignevector[p][sortindex[1]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u1kdxdy_ket u1kdxdy_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(HaldaneHw,0),1)/o/C u2kdxdy_ket u2kdxdy_ket[] = eignevector[p][sortindex[1]] //V_min) lc = W_coefw[2]-0.5*tt*sigma else lc =V_min endif if (W_coefw[2]+0.5*tt*sigma < V_max) lh = W_coefw[2]+0.5*tt*sigma else lh =V_max endif ModifyImage/W=$grabwinnonew(nameofwave(name)) $nameofwave(name) ctab= {lc,lh,VioletOrangeYellow,0} //PRINT num2str(W_coefw[2]-1.5*sigma)+" "+ num2str(W_coefw[2]+1.5*sigma) //print num2str(W_coefw[2]) //display namehistw end Function gethistgram_npcolor(name) string name string histn = "Hist_"+name make/o/N=1000 $histn wavestats/Q $name variable Vmax = V_max Variable Vmin = V_min variable mdian = median($name) Histogram/B={(mdian-(Vmax-Vmin)/2),(V_max-V_min)/1000,1000} $name,$histn duplicate/o $histn fithttemp //CurveFit/Q gauss, $histn/D=fithttemp variable V_FitError = 0 CurveFit/Q gauss, $histn/D=fithttemp if (V_FitError != 0) endif killwaves fithttemp string W_coef="W_coef" wave W_coefw = $W_coef variable coverratio //(0,1] coverratio = (2*5*W_coefw[3])/(Vmax-Vmin) Histogram/B={(mdian-coverratio*(Vmax-Vmin)/2),coverratio*(V_max-V_min)/200,200} $name,$histn //CurveFit/Q gauss, $histn/D//=fithttemp V_FitError = 0 CurveFit/Q gauss, $histn/D//=fithttemp if (V_FitError != 0) endif end Function/s grabwinnonew(name) string name string fulllist = WinList("*", ";","WIN:1") string nn,waveong,cmdn,out out = "" variable i for(i=0; iV_min) lc = W_coefw[2]-0.5*tt*sigma else lc =V_min endif if (W_coefw[2]+0.5*tt*sigma < V_max) lh = W_coefw[2]+0.5*tt*sigma else lh =V_max endif ModifyImage/W=$grabwinnonew(nameofwave(name)) $nameofwave(name) ctab= {lc,lh,VioletOrangeYellow,1} //PRINT num2str(W_coefw[2]-1.5*sigma)+" "+ num2str(W_coefw[2]+1.5*sigma) //print num2str(W_coefw[2]) //display namehistw end Function/s grabwin(name) string name string fulllist = WinList("*", ";","WIN:1") string nn,waveong,cmdn,out out = "" variable i for(i=0; i= 0 && round(qq) <= Ny-1 ) //** condition of Q in range pps = i break endif i+=1 while(i < Nx) return pps end //***************************************************** //## 03 //** Determine the Last P for assignment loop Function findendppf(mat,angle,Zn,addY) String mat variable angle variable Zn variable addY //** Auto ordering layer index make/o/N=3 ordernum ={0,1,2} ordernum=abs(ordernum-zn) wavestats/Q ordernum ordernum={0,1,2} DeletePoints V_minloc,1, ordernum wavestats/Q ordernum variable xn = ordernum[V_minloc] variable yn = ordernum[V_maxloc] killwaves ordernum //print num2str(xn)+" "+num2str(yn) //** Initial layer dimsize wave matw = $mat variable Nx,Ny,Nz Nx = dimsize(matw,xn) Ny = dimsize(matw,yn) Nz = dimsize(matw,zn) //** Search for the End P variable pp, qq,ppe variable i ppe=nan i=Nx-1 //## Search from the largest P, calculate when Q is in range, that position is ppe do qq = tan(angle*pi/180)*(i-round(Nx/2))+round(Ny/2)+addY //[Formula of Rotated HLine] //print qq if(round(qq) <= Ny-1 && round(qq) >= 0) //** condition of Q in range ppe = i break endif i-=1 while(i >= 0) return ppe end //***************************************************** //## 04 //** Find the addY Range for rotated H-Linecut Function/Wave findrangeforangle_LHf(mat,angle,Zn) String mat variable angle variable Zn //** Auto ordering layer index make/o/N=3 ordernum ={0,1,2} ordernum=abs(ordernum-zn) wavestats/Q ordernum ordernum={0,1,2} DeletePoints V_minloc,1, ordernum wavestats/Q ordernum variable xn = ordernum[V_minloc] variable yn = ordernum[V_maxloc] killwaves ordernum //** Search for high limit variable addY, addYmax, addYmin,body addY =0 do body = mod(Round(findendppf(mat,angle,Zn,addY)),1) if (body != 0) addYmax = addY-1 break endif addY+=1 while(addY < 3*dimsize($mat,yn)) //** Search for low limit addY =0 do body = mod(Round(findendppf(mat,angle,Zn,addY)),1) if (body != 0) addYmin = addY+1 break endif addY-=1 while(addY > -3*dimsize($mat,yn)) //** Return the limit wave make/N=2/o Hcut_Vrange //print addYmin //print addYmax Hcut_Vrange={addYmin,addYmax} return Hcut_Vrange end //////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////// //** INTERACTING FUNCTIONAL: Extract rotated Vertical Linecut //////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////////////////// //%%%%%%%%%%%%** get rotated V-linecut, addX here is counted from center. //%%%%%%%%%%%%** V-linecut can move in X direction, addX is DP(Int) //////////////////////////////////////////////////////////////////////////////////////////////////////////// //***************************************************** //## 01 //** Main Function Function anglelinecutVf(mat,angle,Zn,addX,normornot,smornot) String mat variable angle variable Zn variable addX variable normornot variable smornot //** Auto ordering layer index make/o/N=3 ordernum ={0,1,2} ordernum=abs(ordernum-zn) wavestats/Q ordernum ordernum={0,1,2} DeletePoints V_minloc,1, ordernum wavestats/Q ordernum variable xn = ordernum[V_minloc] variable yn = ordernum[V_maxloc] killwaves ordernum //print num2str(xn)+" "+num2str(yn) //** Initial layer dimsize wave matw = $mat variable Nx,Ny,Nz Nx = dimsize(matw,xn) Ny = dimsize(matw,yn) Nz = dimsize(matw,zn) variable pp, qq //** Determine Limit for Assignment Loop //## Play with vertical Line //## Logistically, we swap P and Q //## that mean, set Q as independent variable //## and check P the dependent variable //** Determine the First and Last Q for value assignment loop variable qqs,qqe,pps,ppe,xxs,xxe,yys,yye //pp = tan(angle*pi/180)*(qq-round(Ny/2))+round(Nx/2)+addX qqs=findstartqqf(mat,angle,Zn,addX) //First Q qqe=findendqqf(mat,angle,Zn,addX) //Last Q //** Calculate P and X, Y yys = dimoffset(matw,yn)+qqs*dimdelta(matw,yn) yye = dimoffset(matw,yn)+qqe*dimdelta(matw,yn) pps = tan(-angle*pi/180)*(qqs-round(Ny/2))+round(Nx/2)+addX ppe = tan(-angle*pi/180)*(qqe-round(Ny/2))+round(Nx/2)+addX xxs = dimoffset(matw,xn)+pps*dimdelta(matw,xn) xxe = dimoffset(matw,xn)+ppe*dimdelta(matw,xn) //** Calculate linecut scale variable len = sqrt((xxe-xxs)^2+(yye-yys)^2) string linecutV = "LV_"+mat make/N=((qqe-qqs+1),Nz)/o $linecutV wave linecutVw = $linecutV setscale/p y,dimoffset(matw,zn),dimdelta(matw,zn),"",linecutVw setscale/i x,yys,yys+len,"",linecutVw //** Extract value from 3D matrix //** The value assignment loop is runing along Q variable i i=0 qq=qqs do //qq= round(tan(-angle*pi/180)*(pp-round(Nx/2))+round(Ny/2)+addY) pp = round(tan(-angle*pi/180)*(qq-round(Ny/2))+round(Nx/2)+addX) //[Formula of Rotated VLine] //linecutHw[][i] = matw[pp][qq][p] if(zn == 0) linecutVw[i][] = matw[q][pp][qq] endif if(zn == 1) linecutVw[i][] = matw[pp][q][qq] endif if(zn == 2) linecutVw[i][] = matw[pp][qq][q] endif qq+=1 i+=1 while(qq < qqe+1) //di(linecutHw) //** Shift the zero for FFT linecut string slicename = "Zslice_"+mat wave slicenamew = $slicename variable kk = (xxs-xxe)/(yys-yye) variable bb = xxs-kk*yys variable zx,zy zy = -(bb*kk)/(kk^2+1) zx = bb/(kk^2+1) variable zp,zq zp = round((zx-dimoffset(slicenamew,0))/dimdelta(slicenamew,0)) zq = round((zy-dimoffset(slicenamew,1))/dimdelta(slicenamew,1)) variable dimoff = dimoffset(linecutVw,0) variable dimdel = dimdelta(linecutVw,0) variable correcteddimoffset = -(zq-qqs)*dimdel setscale/p x,correcteddimoffset,dimdel,"",linecutVw //** Shift the zero for FFT linecut //string slicename = "Zslice_"+mat //wave slicenamew = $slicename //variable kk = (xxs-xxe)/(yys-yye) //variable bb = xxs-kk*yys //variable zx,zy //zy = -(bb*kk)/(kk^2+1) //zx = bb/(kk^2+1) //variable zp,zq //zp = round((zx-dimoffset(slicenamew,0))/dimdelta(slicenamew,0)) //zq = round((zy-dimoffset(slicenamew,1))/dimdelta(slicenamew,1)) //correcteddimoffset + dimdelta(linecutHw,0)*zp = 0 //variable correcteddimoffset = -dimdelta(slicenamew,1)*zq //setscale/p x,correcteddimoffset,dimdelta(slicenamew,1),"",linecutVw end //***************************************************** //## 02 //** Determine the First Q for assignment loop Function findstartqqf(mat,angle,Zn,addX) String mat variable angle variable Zn variable addX //** Auto ordering layer index make/o/N=3 ordernum ={0,1,2} ordernum=abs(ordernum-zn) wavestats/Q ordernum ordernum={0,1,2} DeletePoints V_minloc,1, ordernum wavestats/Q ordernum variable xn = ordernum[V_minloc] variable yn = ordernum[V_maxloc] killwaves ordernum //print num2str(xn)+" "+num2str(yn) //** Initial layer dimsize wave matw = $mat variable Nx,Ny,Nz Nx = dimsize(matw,xn) Ny = dimsize(matw,yn) Nz = dimsize(matw,zn) //** Search for the Start Q variable pp, qq,qqs variable i qqs=nan i=0 //## Search from the smallest Q [i.e. 0], calculate when P is in range, that position is pps do pp = tan(-angle*pi/180)*(i-round(Ny/2))+round(Nx/2)+addX //[Formula of Rotated VLine] if(round(pp) >= 0 && round(pp) <= Nx-1 ) //Condition for P in range qqs = i break endif i+=1 while(i < Ny) return qqs end //***************************************************** //## 03 //** Determine the Last Q for assignment loop Function findendqqf(mat,angle,Zn,addX) String mat variable angle variable Zn variable addX //** Auto ordering layer index make/o/N=3 ordernum ={0,1,2} ordernum=abs(ordernum-zn) wavestats/Q ordernum ordernum={0,1,2} DeletePoints V_minloc,1, ordernum wavestats/Q ordernum variable xn = ordernum[V_minloc] variable yn = ordernum[V_maxloc] killwaves ordernum //print num2str(xn)+" "+num2str(yn) //** Initial layer dimsize wave matw = $mat variable Nx,Ny,Nz Nx = dimsize(matw,xn) Ny = dimsize(matw,yn) Nz = dimsize(matw,zn) //** Search for the End Q variable pp, qq,qqe variable i qqe=nan i=Nx-1 //## Search from the largest Q, calculate when P is in range, that position is pps do pp = tan(-angle*pi/180)*(i-round(Ny/2))+round(Nx/2)+addX //[Formula of Rotated VLine] //print qq if(round(pp) >= 0 && round(pp) <= Nx-1) //Condition for P in range qqe = i break endif i-=1 while(i >= 0) return qqe end //***************************************************** //## 04 //** Find the addX Range for rotated V-Linecut Function/Wave findrangeforangle_LVf(mat,angle,Zn) String mat variable angle variable Zn //** Auto ordering layer index make/o/N=3 ordernum ={0,1,2} ordernum=abs(ordernum-zn) wavestats/Q ordernum ordernum={0,1,2} DeletePoints V_minloc,1, ordernum wavestats/Q ordernum variable xn = ordernum[V_minloc] variable yn = ordernum[V_maxloc] killwaves ordernum //** High Limit variable addX, addXmax, addXmin,body addX =0 do body = mod(Round(findendqqf(mat,angle,Zn,addX)),1) if (body != 0) addXmax = addX-1 break endif addX+=1 while(addX < 3*dimsize($mat,xn)) //** Low Limit addX =0 do body = mod(Round(findendppf(mat,angle,Zn,addX)),1) if (body != 0) addXmin = addX+1 break endif addX-=1 while(addX > -3*dimsize($mat,xn)) //** Return Limit make/N=2/o Vcut_Vrange //print addYmin //print addYmax Vcut_Vrange={addXmin,addXmax} return Vcut_Vrange end //////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////// //** Auxiliary Functions //////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////// Function sumoned(dest,orin,zn) string dest //Output name string orin //Input 3D matrix variable zn //Index of your interesting dimension that keep in the end. //** Auto ordering layer index make/o/N=3 ordernum ={0,1,2} ordernum=abs(ordernum-zn) wavestats/Q ordernum ordernum={0,1,2} DeletePoints V_minloc,1, ordernum wavestats/Q ordernum variable xn = ordernum[V_minloc] variable yn = ordernum[V_maxloc] killwaves ordernum wave orinw = $orin make/N=(dimsize(orinw,zn))/o $dest wave destw = $dest destw=0 variable i,j,k,addk,countnum=0 variable is,ie,js,je variable/G FFTmarque_3dplot// 1 for whole, 2 for in marque if (FFTmarque_3dplot == 1) is = 0 ie = dimsize(orinw,xn) js = 0 je = dimsize(orinw,yn) endif if (FFTmarque_3dplot == 2) getmarquee/W=$winname(0,1) left, bottom variable p1,p2,q1,q2 is = round((v_left-dimoffset($tpw(),0))/dimdelta($tpw(),0)) ie = round((v_right-dimoffset($tpw(),0))/dimdelta($tpw(),0))+1 js = round((v_bottom-dimoffset($tpw(),1))/dimdelta($tpw(),1)) je = round((v_top-dimoffset($tpw(),1))/dimdelta($tpw(),1))+1 //print is, ie ,js, je //print tpw() endif //if (mod(Round(mean(orinw)),1)!=0) i=is//round(dimsize(orinw,xn)/6) do j=js//round(dimsize(orinw,yn)/6) do if(zn== 0) destw[]+=orinw[p][i][j] endif if(zn== 1) destw[]+=orinw[i][p][j] endif if(zn== 2) //destw[]+=orinw[i][j][p] k=0 addk=0 do addk+=orinw[i][j][k] k+=1 while(k jumpjudge) //print i InsertPoints dimsize(jp_p1d,0),1, jp_p1d //print dimsize(jp_p1d,0)-1 jp_p1d[dimsize(jp_p1d,0)-1] = (i+1)*sign(deltawavew[i]) endif i+=1 while (i