//Written by Andrew Nelson, anz@ansto.gov.au //last updated Jun05 //uses the Parratt formalism to calculate the reflectivity. //Why is it Motofit, well, Iszy's parrot is called Moto, so there you go. Menu "Reflectivity" "Reflectivity Fit Panel",plotCalcref() "Load experimental data",plotexpref() "Make global coefficients for multiple contrasts",globalsmaker() End Proc plotCalcref(num,qmin,qmax,res,SLDplotnum) Variable num=228,qmin=0.005,qmax=0.3,res=5,SLDplotnum=500 Prompt num "Enter number of data points for model: " Prompt qmin "Enter minimum q-value (A^-1) for model: " Prompt qmax "Enter maximum q-value (A^-1) for model: " Prompt res "Enter %resolution (dq/q): " Prompt SLDplotnum "How many points do you want in the SLD plot" Variable/G plotyp,SLDpts plotyp=1 SLDpts=SLDplotnum Variable logg if(plotyp==1) logg=0 endif if(plotyp==2) logg=1 endif if(plotyp==3) logg=0 endif make/o/d/n=(num) xwave_Cref,ywave_Cref make/o/d/n=(SLDpts) zed,sld //This section pulls up a graph with the reflectivity in it, and a table containing the reflectivity parameters xwave_Cref = (x+1)*((qmax-qmin)/num) make/o/d/n=10 coef_Cref make/o/d/n=1 resolution coef_Cref[0]=1 coef_Cref[1]=1 coef_Cref[2]=0 coef_Cref[3]=2.07 coef_Cref[4]=1e-5 coef_Cref[5]=4 coef_Cref[6]=25 coef_Cref[7]=3.47 coef_Cref[8]=0 coef_Cref[9]=4 make/o/t par_res= {"%resolution (dq/q)"} resolution[0]=res make/o/t/n=13 parameters_Cref parameters_Cref[0]="Numlayers" parameters_Cref[1]="scale" parameters_Cref[2]="SLDtop" parameters_Cref[3]="SLDbase" parameters_Cref[4]="bkg" parameters_Cref[5]="sigma_base" parameters_Cref[6]="thick1" parameters_Cref[7]="SLD1" parameters_Cref[8]="solv1" parameters_Cref[9]="rough1" parameters_Cref[10]="thick2" parameters_Cref[11]="SLD2" parameters_Cref[12]="solv2" parameters_Cref[13]="rough2" //although they don't have to be displayed coz of the panel, I am so that they can be pasted Edit/K=1 parameters_Cref,coef_Cref,par_res,resolution ywave_Cref:= Motofit(coef_Cref,xwave_Cref) Display/K=1 ywave_Cref vs xwave_Cref; AppendToGraph/L=SLDprofile/B=z sld vs zed ModifyGraph standoff(SLDprofile)=0,standoff(z)=0,axisEnab(SLDprofile)={0.6,1} ModifyGraph axisEnab(z)={0.6,1},freePos(SLDprofile)={0,z} ModifyGraph freePos(z)={0,SLDprofile} ModifyGraph log(bottom)=0,mode=0 ModifyGraph log(left)=(logg) Label bottom "q (A\\S-1\\M)" Label left "Reflectivity" //This section pulls up a graph of the (real SLD profile). // It automatically updates whenever you change the fit parameters SLD:= SLDplot(coef_Cref,zed) Label SLDprofile "\Z08\\f02\\F'Symbol'r\\F'Arial' \\f00/10\\S-6\\M ‰\\S-2 " ModifyGraph lblPos(SLDprofile)=42; Modifygraph fSize(SLDprofile)=8 Label z "\Z08z /‰" ModifyGraph lblPos(z)=30; Modifygraph fSize(z)=8 Modifygraph rgb(sld)=(0,0,52224) //this puts all the parameters in a nice convenient panel Reflectivity_panel() End Function Motofit(w,x) :Fitfunc Wave w variable x Wave resolution Variable dq,reflectivity,res Variable/G plotyp NVAR/Z resoln //This function calls calcreflectivity, and integrates over a gaussian beam profile. // Knowledge of dq/q is required. //Wouldn't be difficult at all to modify users different resolution function //This is where the user should set up their own SLD profile If(waveexists(resolution)) res=resolution[1] else res=resoln endif duplicate/o w call Wave call call[4]=0 dq=x*(res/100) reflectivity=calcreflectivity(call,x) if(dq>0) reflectivity+=0.135*calcreflectivity(call,x-dq) reflectivity+=0.135*calcreflectivity(call,x+dq) reflectivity+=0.325*calcreflectivity(call,x-(dq*0.75)) reflectivity+=0.325*calcreflectivity(call,x+(dq*0.75)) reflectivity+=0.605*calcreflectivity(call,x-(dq/2)) reflectivity+=0.605*calcreflectivity(call,x+(dq/2)) reflectivity+=0.88*calcreflectivity(call,x-(dq/4)) reflectivity+=0.88*calcreflectivity(call,x+(dq/4)) reflectivity/=4.89 endif reflectivity+=abs(w[4]) if (plotyp==2) reflectivity=reflectivity endif if(plotyp==3) reflectivity=reflectivity*(x^4) endif if(plotyp==1) reflectivity=log(reflectivity) endif Killwaves/Z kzn,rn,rrn,call return reflectivity End Function Calcreflectivity(w,x) //:fitfunc Wave w Variable x NVAR/Z plotyp Variable reflectivity,ii,nlayers,inter,SLD,qq,scale,bkg,subrough Variable/C super,sub,arg,cinter //this is the kernel that calculates the reflectivity. It works, so don't alter it. //number of layers,SUPERphaseSLD,SUBphaseSLD,Q value //subsequent layers have 4 parameters each: thickness, SLD, solvent penetration and roughness //if you increase the number of layers you have to put extra parameters in. //you should be able to remember the order in which they go. //enter the SLD as SLD*1e6 (it's easier to type). //Layer 1 is always closest to the SUPERPHASE (e.g. air). increasing layers go down //towards the subphase. This may be confusing if you switch between air-solid and solid-liquid //I will write some functions to create exotic SLD profiles if required. nlayers=w[0] scale=w[1] super=w[2]*1e-6 sub=w[3]*1e-6 bkg=abs(w[4]) subrough=w[5] qq=x //for definitions of these see Parratt paper Make/o/d/C/n=(nlayers+2) kzn Make/o/d/C/n=(nlayers+2) rn Make/o/d/C/n=(nlayers+2) RRN //workout the wavevector in the incident medium/superphase inter=cabs(sqrt((qq/2)^2)) kzn[0]=cmplx(inter,0) //workout the wavevector in the subphase kzn[nlayers+1]=sqrt(kzn[0]^2-4*Pi*(sub-super)) //workout the wavevector in each of the layers ii=1 if(ii-1) //reflectivity=abs(Ro)^2 reflectivity=magsqr(RRN[0]) reflectivity*=scale reflectivity+=bkg reflectivity=(reflectivity) return reflectivity End Function SLDplot(w,z) Wave w Variable z NVAR/Z SLDpts Wave zed variable nlayers,SLD1,SLD2,zstart,zend,ii,temp,zinc,summ Variable deltarho,zi,dindex,sigma,thick,dist,rhosolv // //This function calculates the SLD profile. It updates whenever the fitparameters update. // nlayers=w[0] rhosolv=w[3] //setup the start and finish points of the SLD profile if (nlayers==0) zstart=-5-4*abs(w[5]) else zstart=-5-4*abs(w[9]) endif ii=1 temp=0 if (nlayers==0) zend=5+4*abs(w[5]) else do temp+=abs(w[4*ii+2]) ii+=1 while(iinlayers) //kill the parameters ii=oldlayers do index=num2str(ii) controlname="thick"+index killcontrol $controlname controlname="SLD"+index killcontrol $controlname controlname="solv"+index killcontrol $controlname controlname="rough"+index killcontrol $controlname jj=0 parindex=ii*4+2 do //kill the holdwaves index=num2str(parindex) controlname="h"+index killcontrol $controlname parindex+=1 jj+=1 while(jj<4) ii-=1 while (ii>nlayers) elseif(nlayers>oldlayers) //make the new waves Variable ypos //position of the 1st layer String partitle,holdname ii=oldlayers do ii+=1 //ii is the layer number parindex=4*ii+2 //parindex is the parameter number ypos=220+24*(ii-1) //setup the new layer positions partitle=num2istr(parindex)+". thick"+num2istr(ii) //par title is the name of the title controlname="thick"+num2istr(ii) //controlname is the name of the control SetVariable $controlname,pos={31,ypos},size={110,16},title=partitle SetVariable $controlname,value=coef_Cref[parindex] holdname="h"+num2istr(parindex) CheckBox $holdname,pos={144,ypos},size={16,14},title="",value= 0 parindex+=1 partitle=num2istr(parindex)+". SLD"+num2istr(ii) controlname="SLD"+num2str(ii) SetVariable $controlname,pos={172,ypos},size={110,16},title=partitle SetVariable $controlname,value=coef_Cref[parindex] SetVariable $controlname,limits={-Inf,Inf,0.01} holdname="h"+num2istr(parindex) CheckBox $holdname,pos={286,ypos},size={16,14},title="",value= 0 parindex+=1 partitle=num2istr(parindex)+". solv"+num2istr(ii) controlname="solv"+num2str(ii) SetVariable $controlname,pos={310,ypos},size={110,16},title=partitle SetVariable $controlname,value=coef_Cref[parindex] SetVariable $controlname,limits={0,100,1} holdname="h"+num2istr(parindex) CheckBox $holdname,pos={424,ypos},size={16,14},title="",value= 0 parindex+=1 partitle=num2istr(parindex)+". rough"+num2istr(ii) controlname="rough"+num2str(ii) SetVariable $controlname,pos={445,ypos},size={110,16},title=partitle SetVariable $controlname,value=coef_Cref[parindex] SetVariable $controlname,limits={0,100,0.5} holdname="h"+num2istr(parindex) CheckBox $holdname,pos={559,ypos},size={16,14},title="",value= 0 while(ii0) abort "The y data does not have the same number of points as the x data." endif if(waveexists($e)==1) elength=Dimsize($e,0) if(abs(xlength-elength)>0) abort "the error wave does not have the same number of points as the x data" endif endif //do we want to fit with the cursors or not Variable cursors controlinfo fitcursors cursors=V_value String fitdestination="fit_"+y,fitx="fit_"+x Make/o/d/n=(ylength) $fitdestination Duplicate/o $x $fitx //if the cursors aren't on the graph do the fit if(cursors==0) if(waveexists($e)) FuncFit/H=holdstring/N Motofit coef_Cref $y /X=$x /D=$fitdestination /W=$e /I=1 else FuncFit/H=holdstring/N Motofit coef_Cref $y /X=$x /D=$fitdestination endif //here we want to do the fit with the cursors, but we first check that the cursors are on the graph. else if (WaveExists(CsrWaveRef(A)) %& WaveExists(CsrWaveRef(B))) if (CmpStr(CsrWave(A),CsrWave(B)) != 0) abort "The cursors are not on the same wave. Please move them so that they are." endif else abort "The cursors must be placed on the top graph. Select Show Info from the Graph menu for access to the cursors." endif if(waveexists($e)) FuncFit/H=holdstring/N Motofit coef_Cref $y[pcsr(A),pcsr(B)] /X=$x /D=$fitdestination /W=$e /I=1 else FuncFit/H=holdstring/N Motofit coef_Cref $y[pcsr(A),pcsr(B)] /X=$x /D=$fitdestination endif //this cuts off the start and finish of the fit waves if you use the cursor Variable start=pcsr(a),finish=pcsr(b) Deletepoints 0,start, $fitdestination,$fitx Deletepoints (finish-start+1),(ylength-finish-1),$fitdestination,$fitx endif //make the SLD curves and fit curves the same colour variable rr = abs(trunc(enoise(65535))) variable gg = abs(trunc(enoise(65535))) variable bb = abs(trunc(enoise(65535))) //if the trace isn't on the graph add it. If it is, don't do anything. String traceexists=TraceNameList("",";",1) if(Strsearch(traceexists,fitdestination,0)==-1) AppendToGraph $fitdestination vs $fitX ModifyGraph rgb ($fitdestination) =(rr,gg,bb) Modifygraph lsize($fitdestination)=2 endif //figure out what the SLD and Zed and coef waves should be called. Use a stub of the //data wave. unfortunately this means oyu have to go the route of loading data, otherwise //you don't get the logR ending. y=y[0,strlen(y)-2] Doupdate //append an SLD wave as well String SLDdestination="SLD_"+y,zeddestination="zed_"+y Wave sld,zed Duplicate/O sld,$SLDdestination Duplicate/O zed,$zeddestination Setformula $SLDdestination,"" Setformula $zeddestination,"" if(Strsearch(traceexists,SLDdestination,0)==-1) AppendToGraph/L=SLDprofile/B=z $SLDdestination vs $zeddestination ModifyGraph rgb ($SLDdestination) =(rr,gg,bb) Modifygraph lsize($SLDdestination)=2 endif //make a coefficient wave related to the data //y is the name of the y wave data test="coef_"+y Duplicate/O coef_Cref $test End Function Loaddatas(loaddatas) : ButtonControl String Loaddatas NVAR/Z plotyp //this function loads the data into IGOR. However, you will probably have to append the traces to the graph yourself String cmd cmd= "Plotexpref()" Execute/Z cmd End Function Plottype(ctrlName,popNum,popStr) : PopupMenuControl String ctrlName Variable popNum String popStr //this function is an attempt to change the plots without having to kill and reload the data again. //however it doesn't do the fitted waves yet!!! controlinfo plottype Variable olddatatype,newdatatype,ii,ylength olddatatype=plotyp newdatatype=V_Value NVAR/z plotyp plotyp=newdatatype String datalogR,dataR,dataRQ4,temp,ywave,xwave,ewave dataR="" dataR=Wavelist("*R",";","") //if you are going to R vs Q if(newdatatype==2) ii=0 if(olddatatype==1) do temp=stringfromlist(ii,dataR) if (cmpstr(temp[0,3],"coef")==0) else if(strlen(temp)>0) ylength=strlen(temp) ywave=temp[0,ylength-2] xwave=cleanupname((ywave+"q"),0) ewave=cleanupname((ywave+"E"),0) ywave=cleanupname((ywave+"R"),0) Rename $temp, $ywave Wave w=$ywave Wave q=$xwave Wave e=$ewave e=alog(e) //convert to R vs Q from logR w=alog(w) e*=w e-=w endif endif ii+=1 while(strlen(temp)>0) endif ii=0 if(olddatatype==3) do temp=stringfromlist(ii,dataR) if (cmpstr(temp[0,3],"coef")==0) else if(strlen(temp)>0) ylength=strlen(temp) ywave=temp[0,ylength-4] xwave=cleanupname((ywave+"q"),0) ewave=cleanupname((ywave+"E"),0) ywave=cleanupname((ywave+"R"),0) Rename $temp, $ywave Wave w=$ywave Wave q=$xwave Wave e=$ewave e/=q^4 //convert to R vs Q from RQ^4 w=w/(q^4) endif endif ii+=1 while(strlen(temp)>0) endif endif //if you are going to log R vs Q if(newdatatype==1) ii=0 if(olddatatype==2) do temp=stringfromlist(ii,dataR) if (cmpstr(temp[0,3],"coef")==0) else if(strlen(temp)>0) ylength=strlen(temp) ywave=temp[0,ylength-2] xwave=cleanupname((ywave+"q"),0) ewave=cleanupname((ywave+"E"),0) ywave=cleanupname((ywave+"R"),0) Rename $temp, $ywave Wave w=$ywave Wave q=$xwave Wave e=$ewave e=log((e+w)/w) //convert to logR vs Q from R w=log(w) endif endif ii+=1 while(strlen(temp)>0) endif ii=0 if(olddatatype==3) do temp=stringfromlist(ii,dataR) if (cmpstr(temp[0,3],"coef")==0) else if(strlen(temp)>0) ylength=strlen(temp) ywave=temp[0,ylength-2] xwave=cleanupname((ywave+"q"),0) ewave=cleanupname((ywave+"E"),0) ywave=cleanupname((ywave+"R"),0) Rename $temp, $ywave Wave w=$ywave Wave q=$xwave Wave e=$ewave e/=q^4 //convert RQ4 to R w=w/(q^4) e=log((e+w)/w) //convert R to log R w=log(w) endif endif ii+=1 while(strlen(temp)>0) endif endif //if you are going to RQ4 vs Q if(newdatatype==3) ii=0 if(olddatatype==2) do temp=stringfromlist(ii,dataR) if (cmpstr(temp[0,3],"coef")==0) else if(strlen(temp)>0) ylength=strlen(temp) ywave=temp[0,ylength-2] xwave=cleanupname((ywave+"q"),0) ewave=cleanupname((ywave+"E"),0) ywave=cleanupname((ywave+"R"),0) Rename $temp, $ywave Wave w=$ywave Wave q=$xwave Wave e=$ewave e*=q^4 w=w*q^4 endif endif ii+=1 while(strlen(temp)>0) endif ii=0 if(olddatatype==1) do temp=stringfromlist(ii,dataR) if (cmpstr(temp[0,3],"coef")==0) else if(strlen(temp)>0) ylength=strlen(temp) ywave=temp[0,ylength-2] xwave=cleanupname((ywave+"q"),0) ewave=cleanupname((ywave+"E"),0) ywave=cleanupname((ywave+"R"),0) Rename $temp, $ywave Wave w=$ywave Wave q=$xwave Wave e=$ewave e=alog(e) //convert to R vs Q from logR w=alog(w) e*=w e-=w e*=q^4 w=w*q^4 //convert from R to RQ4 endif endif ii+=1 while(strlen(temp)>0) endif endif if(plotyp==2) Modifygraph log(left)=1 else modifygraph log(left)=0 endif End Window globalsmaker() : Panel PauseUpdate; Silent 1 // building window... NewPanel /W=(93.75,355.25,702,676.25) as "Make global coefficient wave" ModifyPanel cbRGB=(0,52224,52224) SetDrawLayer UserBack DrawText 32,57,"0. Number of layers" DrawText 94,77,"1. scale" DrawText 81,98,"2. SLD top" DrawText 69,118,"3. SLD base" DrawText 104,138,"4. bkg" DrawText 73,159,"5. Sig_base" DrawText 40,206,"6. thick1" DrawText 33,226,"10. thick2" DrawText 33,248,"14. thick3" DrawText 33,269,"18. thick4" DrawText 33,291,"22. thick5" DrawText 130,206,"7. SLD1" DrawText 123,228,"11. SLD2" DrawText 123,251,"15. SLD3" DrawText 123,271,"19. SLD4" DrawText 123,293,"23. SLD5" DrawText 217,205,"8. solv1" DrawText 210,225,"12. solv2" DrawText 210,247,"16. solv3" DrawText 210,271,"20. solv4" DrawText 210,290,"24. solv5" DrawText 301,205,"9. rough1" DrawText 294,225,"13. rough2" DrawText 297,248,"17.rough3" DrawText 294,269,"21. rough4" DrawText 294,290,"25. rough5" DrawText 428,169,"coefficient waves" SetDrawEnv fstyle= 1,textrgb= (65280,16384,16384) DrawText 130,34,"global?" SetDrawEnv fstyle= 1,textrgb= (65280,0,0) DrawText 14,34,"coef number" PopupMenu con1,pos={399,201},size={146,21},title="2nd contrast" PopupMenu con1,mode=4,popvalue="_none_",value= #"WaveList(\"coef*\",\";\",\"\") +\"_none_\"" PopupMenu con2,pos={399,228},size={143,21},title="3rd contrast" PopupMenu con2,mode=27,popvalue="_none_",value= #"WaveList(\"coef*\",\";\",\"\") +\"_none_\"" PopupMenu con0,pos={399,175},size={142,21},title="1st contrast" PopupMenu con0,mode=4,popvalue="_none_",value= #"WaveList(\"coef*\",\";\",\"\")+\"_none_\" " PopupMenu con3,pos={399,254},size={143,21},title="4th contrast" PopupMenu con3,mode=4,popvalue="_none_",value= #"WaveList(\"coef*\",\";\",\"\") +\"_none_\"" PopupMenu con4,pos={399,281},size={143,21},title="5th contrast" PopupMenu con4,mode=4,popvalue="_none_",value= #"WaveList(\"coef*\",\";\",\"\") +\"_none_\"" Button makeglobal,pos={389,31},size={162,76},proc=Makeglobals,title="Make global coefficient wave" CheckBox g1,pos={146,63},size={16,14},title="",value= 0 CheckBox g0,pos={146,43},size={16,14},disable=2,title="",value= 1 CheckBox g2,pos={146,84},size={16,14},title="",value= 0 CheckBox g3,pos={146,104},size={16,14},title="",value= 0 CheckBox g4,pos={146,124},size={16,14},title="",value= 0 CheckBox g5,pos={146,146},size={16,14},title="",value= 0 CheckBox g6,pos={91,191},size={16,14},title="",value= 0 CheckBox g10,pos={91,212},size={16,14},title="",value= 0 CheckBox g14,pos={91,234},size={16,14},title="",value= 0 CheckBox g18,pos={91,256},size={16,14},title="",value= 0 CheckBox g22,pos={91,278},size={16,14},title="",value= 0 CheckBox g23,pos={179,278},size={16,14},title="",value= 0 CheckBox g19,pos={179,257},size={16,14},title="",value= 0 CheckBox g15,pos={179,235},size={16,14},title="",value= 0 CheckBox g11,pos={179,212},size={16,14},title="",value= 0 CheckBox g7,pos={179,191},size={16,14},title="",value= 0 CheckBox g8,pos={264,190},size={16,14},title="",value= 0 CheckBox g12,pos={264,211},size={16,14},title="",value= 0 CheckBox g16,pos={264,234},size={16,14},title="",value= 0 CheckBox g20,pos={264,255},size={16,14},title="",value= 0 CheckBox g24,pos={264,277},size={16,14},title="",value= 0 CheckBox g25,pos={359,277},size={16,14},title="",value= 0 CheckBox g21,pos={359,255},size={16,14},title="",value= 0 CheckBox g17,pos={359,233},size={16,14},title="",value= 0 CheckBox g13,pos={359,211},size={16,14},title="",value= 0 CheckBox g9,pos={359,190},size={16,14},title="",value= 0 PopupMenu numlayers,pos={170,39},size={134,21},title="number of layers" PopupMenu numlayers,mode=1,popvalue="0",value= #"\"0;1;2;3;4;5\"" PopupMenu numcontrasts,pos={389,116},size={153,21},title=" number of contrasts" PopupMenu numcontrasts,mode=1,popvalue="1",value= #"\"1;2;3;4;5\"" EndMacro Function Makeglobals(ctrlName) : ButtonControl String ctrlName Variable contrasts,nlayers,ii,npars,numfixed,totalpars,count,jj,kk controlinfo numcontrasts contrasts=V_value controlinfo numlayers nlayers=V_value-1 //V_Value starts from 1, but the number of layers goes from 0. String test,test1,test2,globalstring npars=4*nlayers+6 globalstring="" numfixed=0 //this section makes sure that the waves are the right length (=npars), and are proper waves ii=0 do test2="" test=num2str(ii) test1=cleanupname(("con"+test),0) controlinfo $test1 test2=S_Value if(waveexists($test2)==0) abort "one of the waves might not a coefficient wave, please check" endif if(abs(Dimsize($test2,0)-npars)>0) abort "one of the coefficient waves is the wrong length for the number of layers" endif ii+=1 while(ii