4Misc_StartP#4Platform@9VersionCheck xHH@Rg(HHdh xHH@Rg(HHdh x HH@Rg(HHdh ^Graph*]WDashSettings#  !4 4 4 4 4 4 whomeCCdMacintosh HD:Users:lingyuan:Desktop:RecentWindowsOAuxiliaryd3dErrors.ihfHaldane ModelQi-Hughes-Zhang ModelQi-Wu-Zhang Model 4Misc_EndP#4XOPState_StartP#4XOPState_EndP#\k AYCM+<!*Kd TXET0RGIaKQi-Hughes-Zhang ModelQi-Hughes-Zhang Model#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 ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// //****************** Qi-Hughes-Zhang Model ************************** ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// //Chiral topological superconductor from the quantum Hall state //Xiao-Liang Qi, Taylor L. Hughes, and Shou-Cheng Zhang //Phys. Rev. B 82, 184516 – Published 10 November 2010 //https://journals.aps.org/prb/abstract/10.1103/PhysRevB.82.184516 ///////////////////////////////////////////////////////////////////// //** Text the form of Hamiltonian //** We find the current form missing the PH0(k)P^-1 = -conj(H0(-k)) Function QHZT() //** Create Parameter waves wT("a1",{{"A*sinX"}}) wT("a2",{{"A*sinY"}}) wT("a3",{{"m+B*(2-cosX-cosY)"}}) wT("a4",{{"-u"}}) wT("a5",{{"-∆"}}) wT("a6",{{"-∆2"}}) //wT("e0t",{{"-µ+(kx^2+ky^2)/(2*m)"}}) //wT("Vxyt",{{"a*ky*kx"}}) //wT("Vxt",{{"γso*kx"}}) //wT("Vyt",{{"γso*ky"}}) //wT("mdd0t",{{"-∆2*kx*ky/k0^2"}}) //wT("mddzt",{{"-∆0"}}) //** Create Pauli sequence wave Wn("ps",{{3,1},{3,2},{3,3},{3,0},{2,2},{1,2}}) //** Do text matrix derivation string ps = "ps" automatrixT("a1;a2;a3;a4;a5;a6",$ps) //alternatively you can also do automatrixT("e0t;Vxyt;Vxt;Vyt;mdd0t;mddzt",Wn("ps",{{3,0,0},{3,3,0},{3,1,2},{0,1,1},{2,0,2},{2,3,2}})) end ///////////////////////////////////////////////////////////////////// //** Correct QHZ BdG Hamiltonian Function/Wave QHZ(A,B,m,mu,ReD,ImD,kx,ky) variable A,B,m,mu,ReD,ImD,kx,ky //** Pointer String for called Functions string ps = "ps_QHZ" string result_C = "result_C" //** Create QHZ Hamiltonian variable pf1,pf2,pf3,pf4,pf5,pf6 pf1 = A*sin(kx) pf2 = A*sin(ky) pf3 = m+B*(2-cos(kx)-cos(ky)) pf4 = -mu //** Select a pairing symmetry variable/G pairmode if (pairmode == 1) // S pf5 = -ReD ImD = 0 pf6 = -ImD endif if (pairmode == 2) // S1+iS2 pf5 = -ReD pf6 = -ImD endif if (pairmode == 3) // Px+iPy pf5 = -ReD*sin(kx) pf6 = -ReD*sin(ky) endif if (pairmode == 4) // S+- pf5 = -ReD*cos(kx)*cos(ky) ImD = 0 pf6 = -ImD endif if (pairmode == 5) // Dx2-y2 pf5 = -ReD*(cos(kx)-cos(ky)) ImD = 0 pf6 = -ImD endif if (pairmode == 6) // Dxy pf5 = -ReD*(sin(kx)*sin(ky)) ImD = 0 pf6 = -ImD endif if (pairmode == 7) // Px pf5 = -ReD*sin(kx) ImD = 0 pf6 = -ImD endif if (pairmode == 8) // Py pf5 = -ReD*sin(ky) ImD = 0 pf6 = -ImD endif wC("pf1_QHZ",{{cmplx(pf1,0)}}) wC("pf2_QHZ",{{cmplx(pf2,0)}}) wC("pf3_QHZ",{{cmplx(pf3,0)}}) wC("pf4_QHZ",{{cmplx(pf4,0)}}) wC("pf5_QHZ",{{cmplx(pf5,0)}}) wC("pf6_QHZ",{{cmplx(pf6,0)}}) Wn("ps_QHZ",{{3,1},{3,2},{3,3},{3,0},{2,2},{1,2}}) automatrixC("pf1_QHZ;pf2_QHZ;pf3_QHZ;pf4_QHZ;pf5_QHZ;pf6_QHZ",$ps) wave/C result_Cw = $result_C duplicate/C/o result_Cw QHZH //This is to make k to -k and apply complex conjugation variable/c e23 = QHZH[2][3] QHZH[2][3] = -conj(e23) variable/c e32 = QHZH[3][2] QHZH[3][2] = -conj(e32) //print e23, QHZH[2][3] killwaves $result_C return QHZH end ///////////////////////////////////////////////////////////////////// //** 3D Spetral Function //QHZSpectralFunction(1,1,0.5,0,1,0,0.1,20) Function QHZSpectralFunction(A,B,m,mu,ReD,ImD,selfimg,numxy) variable A,B,m,mu,ReD,ImD variable selfimg variable numxy //= 100 variable/G zn_cons = 2 //** Make void matrix for spectral function make/o/n=(numxy,numxy,500) modelresult_QHZ setscale/i x,-pi,pi,"",modelresult_QHZ setscale/i y,-pi,pi,"",modelresult_QHZ setscale/i z,-5,5,"",modelresult_QHZ modelresult_QHZ = 0 variable roff = dimoffset(modelresult_QHZ,2) variable rdelta = dimdelta(modelresult_QHZ,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 QHZH = "QHZH" //** Start loop for (kx, ky) i=0 do j=0 do //** Define kx, ky by void matrix shape kx = dimoffset(modelresult_QHZ,0)+i*dimdelta(modelresult_QHZ,0) ky = dimoffset(modelresult_QHZ,1)+j*dimdelta(modelresult_QHZ,1) QHZ(A,B,m,mu,ReD,ImD,kx,ky) wave/C QHZHw = $QHZH //** Solve eigenvalue of QHZ Hamiltonian matrixEigenV QHZHw wave/C n=$W_eigenvalues make/n=(dimsize(QHZHw,0))/o sorteigen sorteigen[0]= real(N[0]) sorteigen[1]= real(N[1]) sorteigen[2]= real(N[2]) sorteigen[3]= real(N[3]) //Make 3D Spectral Function t=0 do modelresult_QHZ[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_QHZ,1)) i+=1 while(i < dimsize(modelresult_QHZ,0) ) //** Call smart 3D displayer duplicate/o modelresult_QHZ Z1_001 killwaves modelresult_QHZ sorteigen QHZHw n execute "d3dsimu(\"Z1_001\",2)" //** Modify the 3D player for this special simulation use (band structure) string aa =winname(0,1)+"#G0"; SetAxis/W=$aa/A/R right 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 tilewindows/WINS=winname(0,1)/R/w=(49,0,83,85)/A=(1,1) tilewindows/WINS=winname(1,1)/R/w=(49,0,83,85)/A=(1,1) //string/G mat3dn_cons //string mat3dn = mat3dn_cons+"_FFT3d" //string mat3dn_consf = mat3dn //variable/G zn_cons //the Energy in which dimension? variable/G z_cons SetVariable setvarz_cons title="Z",size={65,14},value=z_cons,proc=SetVarProc_Cons3dplot2 //SetVariable setvarz_cons limits={dimoffset($mat3dn_consf,zn_cons),(dimoffset($mat3dn_consf,zn_cons)+dimdelta($mat3dn_consf,zn_cons)*(dimsize($mat3dn_consf,zn_cons)-1)),dimdelta($mat3dn_consf,zn_cons)} end ///////////////////////////////////////////////////////////////////// //** Interactive Tune band structure Function ButtonProc_QHZ(ctrlName) : ButtonControl String ctrlName Execute "Cons_QHZcut2c()" end Proc Cons_QHZcut2c(A,B,m,mu,ReD,ImD,selfimg,numxy,mode) variable A = 1 //= A_QHZ variable B = 1//= B_QHZ variable m = 0//= m_QHZ variable mu = 0//= mu_QHZ variable ReD = 0.5//= ReD_QHZ variable ImD = 0//= ImD_QHZ variable selfimg = 0.03 //= selfimg_QHZ variable numxy = 200//= numxy_QHZ variable mode //= pairmode Prompt mu,"μ" prompt ReD,"Re(Δ)" Prompt ImD,"Im(Δ)" Prompt selfimg,"Im(Σ)" prompt numxy,"Number of point in k space" Prompt mode,"Pairing sym",popup,"s;s+is;px+ipy;s±;dx2-y2;dxy;px;py" QHZ_bandtune(A,B,m,mu,ReD,ImD,selfimg,numxy,mode) end Function QHZ_bandtune(A,B,m,mu,ReD,ImD,selfimg,numxy,mode) variable A,B,m,mu,ReD,ImD,selfimg,numxy,mode variable/G pairmode = mode Cons_QHZcut(A,B,m,mu,ReD,ImD,selfimg,numxy) Cons_QHZcut2(A,B,m,mu,ReD,ImD,selfimg,numxy) wave cut_QHZ2=$"cut_QHZ2" wave cut_QHZ=$"cut_QHZ" display/N=QHZmodelinteractive;modifygraph width=700,height=500 Variable/G A_QHZ = A Variable/G B_QHZ = B Variable/G m_QHZ = m Variable/G mu_QHZ = mu Variable/G ReD_QHZ = ReD Variable/G ImD_QHZ = ImD Variable/G selfimg_QHZ = selfimg Variable/G numxy_QHZ = numxy variable/G pairmode = mode SetVariable sett1_QHZ win=QHZmodelinteractive, title="A",size={50,15},value=A_QHZ,limits={-inf,inf,0.1},proc=SetVarProc_QHZ_band SetVariable sett2_QHZ win=QHZmodelinteractive, title="B",size={50,15},value=B_QHZ,limits={-inf,inf,0.1},proc=SetVarProc_QHZ_band,pos={53,1} SetVariable setM_QHZ win=QHZmodelinteractive, title="M",size={50,15},value=m_QHZ,limits={-inf,inf,0.05},proc=SetVarProc_QHZ_band,pos={105,1} SetVariable setphi_QHZ win=QHZmodelinteractive, title="μ",size={50,15},value=mu_QHZ,limits={-inf,inf,0.1},proc=SetVarProc_QHZ_band,pos={158,1} SetVariable setreal_QHZ win=QHZmodelinteractive, title="Re(Δ)",size={70,15},value=ReD_QHZ,limits={-inf,inf,0.1},proc=SetVarProc_QHZ_band,pos={208,1} if(pairmode == 2) SetVariable setimg_QHZ win=QHZmodelinteractive, title="Im(Δ)",size={70,15},value=ImD_QHZ,limits={-inf,inf,0.1},proc=SetVarProc_QHZ_band,pos={278,1} else if(pairmode == 3) SetVariable setimg_QHZ win=QHZmodelinteractive, title="Im(Δ)",size={70,15},value=ReD_QHZ,limits={ReD_QHZ,ReD_QHZ,0},pos={278,1} else SetVariable setimg_QHZ win=QHZmodelinteractive, title="Im(Δ)",size={70,15},value=ImD_QHZ,limits={ImD_QHZ,ImD_QHZ,0},proc=SetVarProc_QHZ_band,pos={278,1} endif endif SetVariable setselfimg_QHZ win=QHZmodelinteractive, title="Im(Σ)",size={70,15},value=selfimg_QHZ,limits={-inf,inf,0.1},proc=SetVarProc_QHZ_band,pos={368,1} SetVariable setnum_QHZ win=QHZmodelinteractive, title="num",size={70,15},value=numxy_QHZ,limits={0,inf,10},proc=SetVarProc_QHZ_band,pos={440,1} Button b3sQHZ win=QHZmodelinteractive, title="3D",size={70,15},proc=ButtonProc_QHZband Button m0BCQHZ win=QHZmodelinteractive, title="Chern(m=0)",size={100,15},proc=ButtonProc_QHZChernm0,pos={81,21},fSize=11 Button m0BCQHZ2 win=QHZmodelinteractive, title="Z(Δ,m=0)",size={100,15},proc=ButtonProc_QHZChernm02,pos={81,40},fSize=11 Button d0BCQHZ win=QHZmodelinteractive, title="Chern(Δ=0)",size={100,15},proc=ButtonProc_QHZChernd0,pos={188,21},fSize=11 Button d0BCQHZ2 win=QHZmodelinteractive, title="Z(m,Δ=0,μ=μ0)",size={100,15},proc=ButtonProc_QHZChernd0l1,pos={188,40},fSize=11 Button d0BCQHZ3 win=QHZmodelinteractive, title="Z(μ,Δ=0,m=m0)",size={100,15},proc=ButtonProc_QHZChernd0l2,pos={188,60},fSize=11 Button Edge win=QHZmodelinteractive, title="Slab",size={50,15},proc=ButtonProc_QHZedge,pos={650,1},fSize=11 PopupMenu gapsym proc=PopMenuProc_QHZgapsym,value="s;s+is;px+ipy;s±;dx2-y2;dxy;px;py",title="Gap. Sym.",pos={516,0},mode=mode Display/HOST=#/W=(0.0,0.05,0.5,1);appendimage cut_QHZ;ModifyGraph mirror=2;ModifyImage cut_QHZ ctab= {*,*,VioletOrangeYellow,0};ModifyGraph zero(left)=8;Label bottom "\\Z16X \t \t\t\t\t\t\t\t\tΓ\t\t\t\t\t\t\t\t\t X" setActiveSubwindow ##;Display/HOST=#/W=(0.5,0.05,1,1);appendimage cut_QHZ2;ModifyGraph mirror=2;ModifyImage cut_QHZ2 ctab= {*,*,VioletOrangeYellow,0};ModifyGraph zero(left)=8;Label bottom "\\Z16M \t \t\t\t\t\t\t\t\tΓ\t\t\t\t\t\t\t\t\t M" setActiveSubwindow ##; end Function PopMenuProc_QHZgapsym(ctrlName,popNum,popStr) : PopupMenuControl String ctrlName Variable popNum String popStr Variable/G pairmode variable/G ImD_QHZ variable/G ReD_QHZ pairmode=popNum if(pairmode == 2) SetVariable setimg_QHZ win=QHZmodelinteractive, title="Im(Δ)",size={70,15},value=ImD_QHZ,limits={-inf,inf,0.1},proc=SetVarProc_QHZ_band,pos={278,1} else if(pairmode == 3) SetVariable setimg_QHZ win=QHZmodelinteractive, title="Im(Δ)",size={70,15},value=ReD_QHZ,limits={ReD_QHZ,ReD_QHZ,0},pos={278,1} else SetVariable setimg_QHZ win=QHZmodelinteractive, title="Im(Δ)",size={70,15},value=ImD_QHZ,limits={ImD_QHZ,ImD_QHZ,0},proc=SetVarProc_QHZ_band,pos={278,1} endif endif Variable/G A_QHZ //= A Variable/G B_QHZ //= B Variable/G m_QHZ //= m Variable/G mu_QHZ //= mu Variable/G selfimg_QHZ //= selfimg Variable/G numxy_QHZ //= numxy variable A = A_QHZ variable B = B_QHZ variable m = m_QHZ variable mu = mu_QHZ variable ReD = ReD_QHZ variable ImD = ImD_QHZ variable selfimg = selfimg_QHZ variable numxy = numxy_QHZ Cons_QHZcut(A,B,m,mu,ReD,ImD,selfimg,numxy) Cons_QHZcut2(A,B,m,mu,ReD,ImD,selfimg,numxy) if (cmpstr(grabwinnonew("phasediagramtempQHZ"),"") == 1) if (pairmode == 2 || pairmode == 6 || pairmode == 8) Solveedgestate_QHZ_cmplx(A_QHZ,B_QHZ,m_QHZ,mu_QHZ,ReD_QHZ,ImD_QHZ,10,300) else Solveedgestate_QHZ(A_QHZ,B_QHZ,m_QHZ,mu_QHZ,ReD_QHZ,ImD_QHZ,10,300) endif endif variable/G LDOSEawr_QHZ if (LDOSEawr_QHZ == 1) effQHZ_LDOSEawr() endif End Function ButtonProc_QHZband(ctrlName) : ButtonControl String ctrlName Variable/G A_QHZ //= A Variable/G B_QHZ //= B Variable/G m_QHZ //= m Variable/G mu_QHZ //= mu Variable/G ReD_QHZ //= ReD Variable/G ImD_QHZ //= ImD Variable/G selfimg_QHZ //= selfimg variable A = A_QHZ variable B = B_QHZ variable m = m_QHZ variable mu = mu_QHZ variable ReD = ReD_QHZ variable ImD = ImD_QHZ variable selfimg = selfimg_QHZ QHZSpectralFunction(A,B,m,mu,ReD,ImD,selfimg,100) QWZSpectralFunction2(A,B,m,mu,selfimg,100) end Function Cons_QHZcut(A,B,m,mu,ReD,ImD,selfimg,numxy) variable A,B,m,mu,ReD,ImD,selfimg,numxy //** Make void matrix for spectral function make/o/n=(numxy,500) cut_QHZ setscale/i x,-pi,pi,"",cut_QHZ setscale/i y,-5,5,"",cut_QHZ cut_QHZ = 0 variable roff = dimoffset(cut_QHZ,1) variable rdelta = dimdelta(cut_QHZ,1) //** Variable for Spectral Function variable i,j,selfreal,t selfreal = 0 variable kx variable ky = 0 //** Pointer String for called Functions string W_eigenvalues="W_eigenvalues" string QHZH = "QHZH" //** Start loop for (kx) i=0 do //** Define kx, ky by void matrix shape kx = dimoffset(cut_QHZ,0)+i*dimdelta(cut_QHZ,0) QHZ(A,B,m,mu,ReD,ImD,kx,ky) wave/C QHZHw = $QHZH //** Solve eigenvalue of QHZ Hamiltonian matrixEigenV QHZHw wave/C n=$W_eigenvalues make/n=(dimsize(QHZHw,0))/o sorteigen sorteigen[0]= real(N[0]) sorteigen[1]= real(N[1]) sorteigen[2]= real(N[2]) sorteigen[3]= real(N[3]) //Make 3D Spectral Function t=0 do cut_QHZ[i][]+=(-1/pi)*selfimg/(selfimg^2+((roff+q*rdelta)-selfreal-sorteigen[t])^2) t+=1 while (t < dimsize(sorteigen,0)) i+=1 while(i < dimsize(cut_QHZ,0) ) killwaves QHZHw sorteigen n //di(cut_QHZ) end Function Cons_QHZcut2(A,B,m,mu,ReD,ImD,selfimg,numxy) variable A,B,m,mu,ReD,ImD,selfimg,numxy //** Make void matrix for spectral function make/o/n=(numxy,500) cut_QHZ2 setscale/i x,-sqrt(2)*pi,sqrt(2)*pi,"",cut_QHZ2 setscale/i y,-5,5,"",cut_QHZ2 cut_QHZ2 = 0 variable roff = dimoffset(cut_QHZ2,1) variable rdelta = dimdelta(cut_QHZ2,1) //** 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 QHZH = "QHZH" //** Start loop for (kx) i=0 do //** Define kx, ky by void matrix shape kx = -pi+i*2*pi/(dimsize(cut_QHZ2,0)-1) ky = kx QHZ(A,B,m,mu,ReD,ImD,kx,ky) wave/C QHZHw = $QHZH //** Solve eigenvalue of QHZ Hamiltonian matrixEigenV QHZHw wave/C n=$W_eigenvalues make/n=(dimsize(QHZHw,0))/o sorteigen sorteigen[0]= real(N[0]) sorteigen[1]= real(N[1]) sorteigen[2]= real(N[2]) sorteigen[3]= real(N[3]) //Make 3D Spectral Function t=0 do cut_QHZ2[i][]+=(-1/pi)*selfimg/(selfimg^2+((roff+q*rdelta)-selfreal-sorteigen[t])^2) t+=1 while (t < dimsize(sorteigen,0)) i+=1 while(i < dimsize(cut_QHZ2,0) ) killwaves QHZHw sorteigen n //di(cut_QHZ) end Function SetVarProc_QHZ_band(ctrlName,varNum,varStr,varName) : SetVariableControl String ctrlName Variable varNum String varStr String varName Variable/G A_QHZ //= A Variable/G B_QHZ //= B Variable/G m_QHZ //= m Variable/G mu_QHZ //= mu Variable/G ReD_QHZ //= ReD Variable/G ImD_QHZ //= ImD Variable/G selfimg_QHZ //= selfimg Variable/G numxy_QHZ //= numxy variable/G pairmode variable A = A_QHZ variable B = B_QHZ variable m = m_QHZ variable mu = mu_QHZ variable ReD = ReD_QHZ variable ImD = ImD_QHZ variable selfimg = selfimg_QHZ variable numxy = numxy_QHZ Cons_QHZcut(A,B,m,mu,ReD,ImD,selfimg,numxy) Cons_QHZcut2(A,B,m,mu,ReD,ImD,selfimg,numxy) //print cmpstr(grabwinnonew("phasediagramtempQHZ"),"") if (cmpstr(grabwinnonew("phasediagramtempQHZ"),"") == 1) if (pairmode == 2 || pairmode == 6 || pairmode == 8) Solveedgestate_QHZ_cmplx(A_QHZ,B_QHZ,m_QHZ,mu_QHZ,ReD_QHZ,ImD_QHZ,10,300) else Solveedgestate_QHZ(A_QHZ,B_QHZ,m_QHZ,mu_QHZ,ReD_QHZ,ImD_QHZ,10,300) endif //Cursor/W=QHZedgemodeshowinteractive/I b phasediagramtempQHZ m, ReD endif variable/G LDOSEawr_QHZ if (LDOSEawr_QHZ == 1) effQHZ_LDOSEawr() endif end ///////////////////////////////////////////////////////////////////// //** Calculate Chern number when m = 0 Function ButtonProc_QHZChernm0(ctrlName) : ButtonControl String ctrlName QHZ_hpBerryCurvature(30) QHZ_hmBerryCurvature(30) wave BC_QHZ_hm1 = $"BC_QHZ_hm1" wave BC_QHZ_hm2 = $"BC_QHZ_hm2" wave BC_QHZ_hp1 = $"BC_QHZ_hp1" wave BC_QHZ_hp2 = $"BC_QHZ_hp2" display/N=QHZeq5BSshowfig;modifygraph width= 350,height=300 Display/HOST=#/W=(0,0,0.5,0.5);appendimage BC_QHZ_hm1;ModifyGraph width={Plan,1,bottom,left},mirror=2;ModifyImage BC_QHZ_hm1 ctab= {*,*,VioletOrangeYellow,0};Label left "h-(#1)" setActiveSubwindow ##;Display/HOST=#/W=(0.5,0,1,0.5);appendimage BC_QHZ_hm2;ModifyGraph width={Plan,1,bottom,left},mirror=2;ModifyImage BC_QHZ_hm2 ctab= {*,*,VioletOrangeYellow,0};Label left "h-(#2)" setActiveSubwindow ##;Display/HOST=#/W=(0,0.5,0.5,1);appendimage BC_QHZ_hp1;ModifyGraph width={Plan,1,bottom,left},mirror=2;ModifyImage BC_QHZ_hp1 ctab= {*,*,VioletOrangeYellow,0};Label left "h+(#1)" setActiveSubwindow ##;Display/HOST=#/W=(0.5,0.5,1,1);appendimage BC_QHZ_hp2;ModifyGraph width={Plan,1,bottom,left},mirror=2;ModifyImage BC_QHZ_hp2 ctab= {*,*,VioletOrangeYellow,0};Label left "h+(#2)" setActiveSubwindow ##; ckfig_child(winname(0,1)) tilewindows/WINS=winname(0,1)/R/A=(1,1)/w=(3,58,50,100) end Function ButtonProc_QHZChernm02(ctrlName) : ButtonControl String ctrlName QHZeq8Chernline() end ///////////////////////////////////////////////////////////////////// //# PRB 82, 184516 eq. (5) calculate the chern number of m = 0 case Function qHZeq5hphm(kx, ky) variable kx variable ky Variable/G A_QHZ //= A Variable/G B_QHZ //= B Variable/G ReD_QHZ //= ReD Variable/G ImD_QHZ Variable/G numxy_QHZ //= numxy variable A = A_QHZ variable B = B_QHZ variable ReD = ReD_QHZ variable ImD = ImD_QHZ variable numxy = numxy_QHZ make/o/N=(2,2)/C QHZ_hp QHZ_hp[0][0] = cmplx(sqrt(ReD^2+ImD^2) + B*(2-cos(kx)-cos(ky)),0) QHZ_hp[0][1] = cmplx(A*sin(kx),-A*sin(ky)) QHZ_hp[1][0] = cmplx(A*sin(kx),A*sin(ky)) QHZ_hp[1][1] = -cmplx(sqrt(ReD^2+ImD^2) + B*(2-cos(kx)-cos(ky)),0) make/o/N=(2,2)/C QHZ_hm QHZ_hm[0][0] = cmplx(-sqrt(ReD^2+ImD^2) + B*(2-cos(kx)-cos(ky)),0) QHZ_hm[0][1] = cmplx(A*sin(kx),-A*sin(ky)) QHZ_hm[1][0] = cmplx(A*sin(kx),A*sin(ky)) QHZ_hm[1][1] = -cmplx(-sqrt(ReD^2+ImD^2) + B*(2-cos(kx)-cos(ky)),0) end Function QHZ_hpBerryCurvature(numxy) 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_QHZ_hp setscale/i x,-pi,pi,"",voidberry_QHZ_hp setscale/i y,-pi,pi,"",voidberry_QHZ_hp //setscale/i z,-2.5,4.5,"",voidberry_QHZ_hp voidberry_QHZ_hp = 0 //** Make void matrix for Berry curvature band #1 (lower) and #2 (higher) make/o/n=(numxy,numxy) BC_QHZ_hp1 setscale/i x,-pi,pi,"",BC_QHZ_hp1 setscale/i y,-pi,pi,"",BC_QHZ_hp1 BC_QHZ_hp1 = 0//cmplx(0,0) make/o/n=(numxy,numxy) BC_QHZ_hp2 setscale/i x,-pi,pi,"",BC_QHZ_hp2 setscale/i y,-pi,pi,"",BC_QHZ_hp2 BC_QHZ_hp2 = 0//cmplx(0,0) //** Pointer String for called Functions string W_eigenvalues="W_eigenvalues" string QHZ_hpH = "QHZ_hp" 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_QHZ_hp,0)+i*dimdelta(voidberry_QHZ_hp,0) ky = dimoffset(voidberry_QHZ_hp,1)+j*dimdelta(voidberry_QHZ_hp,1) qHZeq5hphm(kx, ky) wave/C QHZ_hpHw = $QHZ_hpH //** Solve eigenvalue of QHZ_hp Hamiltonian matrixEigenV/R QHZ_hpHw wave/C n=$W_eigenvalues wave/C eignevector=$M_R_eigenVectors make/n=(dimsize(QHZ_hpHw,0))/o sorteigen sorteigen[0]= real(N[0]) sorteigen[1]= real(N[1]) //**Sort eigenvalue and eigenvector make/o/n=(dimsize(QHZ_hpHw,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(QHZ_hpHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZ_hpHw,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(QHZ_hmHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZ_hmHw,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(QHZ_hpHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZ_hpHw,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(QHZ_hmHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZ_hmHw,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(QHZ_d0pHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZ_d0pHw,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(QHZ_d0mHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZ_d0mHw,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(QHZ_d0pHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZ_d0pHw,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(QHZ_d0mHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZ_d0mHw,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] E3k = sorteigen[2] E4k = sorteigen[3] //** make |u1,k> //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u2k_ket u2k_ket[] = eignevector[p][sortindex[1]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u3k_ket u3k_ket[] = eignevector[p][sortindex[2]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u4k_ket u4k_ket[] = eignevector[p][sortindex[3]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u2k_ket u2k_ket[] = eignevector[p][sortindex[1]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u1kdx_ket u1kdx_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u2kdx_ket u2kdx_ket[] = eignevector[p][sortindex[1]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u1kdy_ket u1kdy_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u2kdy_ket u2kdy_ket[] = eignevector[p][sortindex[1]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u1kdxdy_ket u1kdxdy_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QHZHw,0),1)/o/C u2kdxdy_ket u2kdxdy_ket[] = eignevector[p][sortindex[1]] // 0) if (abs(xx) > abs(yy)) phasediagramtempQHZ[i][j] = 1 endif if (abs(xx) < abs(yy)) phasediagramtempQHZ[i][j] = 2 endif else if (abs(xx) > abs(yy)) phasediagramtempQHZ[i][j] = 3 endif if (abs(xx) < abs(yy)) phasediagramtempQHZ[i][j] = 2 endif endif j+=1 while (j sin(kx) pf2 = A*sin(ky) pf3 = m+B*(2-cos(kx)-cos(ky)) //Tranform from m+B*(kx^2+ky^2); kx^2 --> 1-cos(kx) pf4 = -mu wC("pf1_QWZ",{{cmplx(pf1,0)}}) wC("pf2_QWZ",{{cmplx(pf2,0)}}) wC("pf3_QWZ",{{cmplx(pf3,0)}}) wC("pf4_QWZ",{{cmplx(pf4,0)}}) Wn("ps_QWZ",{{1},{2},{3},{0}}) automatrixC("pf1_QWZ;pf2_QWZ;pf3_QWZ;pf4_QWZ",$ps) wave/C result_Cw = $result_C duplicate/C/o result_Cw QWZH killwaves $result_C $ps return QWZH //** Error message, no stop but Continuing execution String msg Variable err msg=GetErrMessage(GetRTError(0),3); err=GetRTError(1) if (err != 0) //Print "Error in Demo: " + msg //Print "Continuing execution" endif end Function QWZSpectralFunction(A,B,m,mu,selfimg,numxy) variable A,B,m,mu variable selfimg variable numxy //= 100 variable/G zn_cons = 2 //** Make void matrix for spectral function make/o/n=(numxy,numxy,500) modelresult_QWZ setscale/i x,-pi,pi,"",modelresult_QWZ setscale/i y,-pi,pi,"",modelresult_QWZ setscale/i z,-3,3,"",modelresult_QWZ modelresult_QWZ = 0 variable roff = dimoffset(modelresult_QWZ,2) variable rdelta = dimdelta(modelresult_QWZ,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 QWZH = "QWZH" //** Start loop for (kx, ky) i=0 do j=0 do //** Define kx, ky by void matrix shape kx = dimoffset(modelresult_QWZ,0)+i*dimdelta(modelresult_QWZ,0) ky = dimoffset(modelresult_QWZ,1)+j*dimdelta(modelresult_QWZ,1) QWZ(A,B,m,kx,ky,mu) wave/C QWZHw = $QWZH //** Solve eigenvalue of QWZ Hamiltonian matrixEigenV QWZHw wave/C n=$W_eigenvalues make/n=(dimsize(QWZHw,0))/o sorteigen sorteigen[0]= real(N[0]) sorteigen[1]= real(N[1]) //Make 3D Spectral Function t=0 do modelresult_QWZ[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_QWZ,1)) i+=1 while(i < dimsize(modelresult_QWZ,0) ) //** Call smart 3D displayer duplicate/o modelresult_QWZ Z2_001 //killwaves modelresult_QWZ sorteigen QWZHw n execute "d3dsimu(\"Z2_001\",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,RedWhiteBlue,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 Function QWZBerryCurvature(A,B,m,mu,numxy) variable A,B,m,mu 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_QWZ setscale/i x,-pi,pi,"",voidberry_QWZ setscale/i y,-pi,pi,"",voidberry_QWZ //setscale/i z,-2.5,4.5,"",voidberry_QWZ voidberry_QWZ = 0 //** Make void matrix for Berry curvature band #1 (lower) and #2 (higher) make/o/n=(numxy,numxy) BC_QWZ1 setscale/i x,-pi,pi,"",BC_QWZ1 setscale/i y,-pi,pi,"",BC_QWZ1 BC_QWZ1 = 0//cmplx(0,0) make/o/n=(numxy,numxy) BC_QWZ2 setscale/i x,-pi,pi,"",BC_QWZ2 setscale/i y,-pi,pi,"",BC_QWZ2 BC_QWZ2 = 0//cmplx(0,0) //** Pointer String for called Functions string W_eigenvalues="W_eigenvalues" string QWZH = "QWZH" 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_QWZ,0)+i*dimdelta(voidberry_QWZ,0) ky = dimoffset(voidberry_QWZ,1)+j*dimdelta(voidberry_QWZ,1) QWZ(A,B,m,kx,ky,mu) wave/C QWZHw = $QWZH //** Solve eigenvalue of QWZ Hamiltonian matrixEigenV/R QWZHw wave/C n=$W_eigenvalues wave/C eignevector=$M_R_eigenVectors make/n=(dimsize(QWZHw,0))/o sorteigen sorteigen[0]= real(N[0]) sorteigen[1]= real(N[1]) //**Sort eigenvalue and eigenvector make/o/n=(dimsize(QWZHw,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(QWZHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QWZHw,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(QWZHw,0),1)/o/C u1k_ket u1k_ket[] = eignevector[p][sortindex[0]] // //|ket> make/n=(dimsize(QWZHw,0),1)/o/C u2k_ket u2k_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