#pragma rtGlobals=1 // Use modern global access method. #pragma version = 2.7 // 2.7 is for changes that occurred in vector-math 2.0 // 2.6 is the Igor 5 version of 2.4 // with version 2.4, I added the offset to bx in the calculation // This is the new version that deals with multiple detecgtors, and has the beam along the z-axis // The traditional beam line coordinates #include "Elements", version>=1.1 #include "Neutron", version>=1.0 #include "vector-math", version>=2.0 // for the new motor used in September, depth = zi - surface - translator Menu "SCD" "行行 General 行行" "tell about spot at cursor A",InfoAboutSpotOnDetector(NaN,NaN,0.13) "Channel to Wavelength",chan2lambda(NaN) "Wavelength to Channel",lambda2chan(Nan) "-" "行行 Reciprocal Lattice 行行" "find hkl of cursor using reip",print hkl_fromDetectorSpot_recip(NaN,NaN,$"") "angle between G vectors at cursors", print angleBetweenGofCursors() "angle between two vectors", angleBetweenHKLs(NaN,NaN,NaN,NaN,NaN,NaN) "detector xy from hkl and recip",printDetectorXY_of_hkl($"",NaN,NaN,NaN) "add predicted hkl spots",addPredictedHKL($"",NaN,NaN) "remove predicted hkl spots", SetDrawLayer /K UserFront // clear the drawing layer "wavelength of hkl from recip",hkl2lambda($"",NaN,NaN,NaN) "hkl pointing at source",hkl_BeamAxis($"") "calc recip lattice from 2 cursors", spots2recip(NaN,NaN,NaN,NaN) "calc recip lattice from 2 cursors xx", xxspots2recip("",NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN) "-" "行行 Displaying 行行" "Show a Slice Step Graph", ShowSliceGraph("") "Run the Slice as a Movie", RunAsMovie() "Apply-Remove gamma Correction",ApplyGamma(NaN) "Make Movie of Top Slice Graph", makeMovieOfTopSliceGraph(NaN) "show SCD geometry",Panel_SCD_geometry() "-" "行行 Sorting and Reading In 行行" "Sort Out 1 Depth", SortDataOneDepth(NaN,NaN,"","",0) "Sum Up Over Wavelength", sumUpWavelength("","",NaN) "Make Sliced Sums at Multiple Depths", MakeSlicedSumsAtManyDepths("","",NaN,NaN,NaN,NaN) "Set One Geometry Parameters",DetectorGeometryWindow() "ReadRangeSCD", ReadRangeSCD("","",nan,"","") "Set Translator Distances", setTranslatorDistances("","",NaN,NaN) "-" "Initialize SCD package",SCDinitPackage(1) End Proc GraphImageStyle() : GraphStyle PauseUpdate; Silent 1 // modifying window... ModifyGraph/Z gfSize=14,width={Aspect,1}, tick=2, minor=1, lowTrip=0.001,standoff=0, mirror=1 ModifyGraph/Z axOffset(left)=-2.3,axOffset(bottom)=-1.11111 ModifyImage $StringFromList(0,ImageNameList("",";")) ctab= {*,*,Terrain,1} // ModifyImage $StringFromList(0,ImageNameList("",";")) ctab= {*,*,Rainbow,0} EndMacro // ==================================================================== // ==================================================================== // ========================= Start of Reciprocal Lattice ======================= // ==================================================================== Function/T hkl_fromDetectorSpot_recip(detX,detY,recip) // computes hkl of a pixel using recip Variable detX,detY Wave/D recip if (!WaveExists(recip) || DimSize(recip,0)!=3 || DimSize(recip,1)!=3 || DimSize(recip,2)!=0) String recipLattice, mat3x3List=WaveList_NxN(3,3) recipLattice = SelectString(exists("recip")==1,"","recip") Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List DoPrompt "reciprocal lattice",recipLattice if (V_flag) return "" endif Wave/D recip = $recipLattice endif detX = numtype(detX) ? hcsr(A) : detX detY = numtype(detY) ? vcsr(A) : detY String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";")))) Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta") Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd") Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx") Variable xp, zp, yp=detY // y is vertical pixelXZfromDetector(deta,detd,detX,xp,zp) // calculate pixel coordinates xp,zp Variable lambda = getLambdaFromPlot("") // get lambda from a step slice plot Variable k_len = 2*PI/lambda Make/N=3/O kf, ki, Gvec ki={0, 0, k_len} // kf={xp,yp,zp} kf={xp-bx,yp,zp} normalize(kf) kf *= k_len // kf and ki need to be k_len long Gvec = kf - ki // recip x H = G --> recipInv x recip x H = recipInv x G // hkl = recipInv x Gvec, so find inverse of recip Make/N=(3,3)/O/D recip_temp_,recipInv recipInv = (p==q) recip_temp_ = recip MatrixLinearSolve /O recip_temp_ recipInv // recipInv now contains the inverse of hklN KillWaves/Z recip_temp_ MatrixMultiply recipInv, Gvec // result is in M_product Wave hkl = M_product String out sprintf out, "%.0f %.0f %.0f",hkl[0],hkl[1],hkl[2] KillWaves/Z ki,kf,Gvec, W_IPIV, recipInv, M_product return out End Function hkl2lambda(recip,h,k,l) // find the wavelength for an hkl using the reciprocal lattice Wave/D recip // 3x3 matrix with reciprocal lattice Variable h,k,l Variable printIt=0 if (numtype(h) || numtype(k) || numtype(l) || !WaveExists(recip)) h = numtype(h) ? 0 : h k = numtype(k) ? 2 : k l = numtype(l) ? 2 : l Prompt h, "h" Prompt k, "k" Prompt l, "l" String recipLattice, mat3x3List=WaveList_NxN(3,3) recipLattice = SelectString(exists("recip")==1,"","recip") Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List DoPrompt "hkl",recipLattice,h,k,l if (V_flag) return NaN endif Wave/D recip = $recipLattice printIt = 1 endif if (numtype(h) || numtype(k) || numtype(l) || !WaveExists(recip)) return NaN endif Make/N=3/D/O hkl__, kihat={0,0,1} hkl__ = {h,k,l} MatrixMultiply recip, hkl__ // result is in M_product KillWaves/Z Ghat Rename M_product Ghat // G vector (from hkl__) Redimension/N=3 Ghat // needed because dimSize(Ghat,1)==1, not 0 Variable G = normalize(Ghat) // d = 2/G // cos(alpha) = -MatrixDot(Ghat,kihat) // theta = /2-alpha // sin(theta) = sin(/2-alpha) = cos(alpha), but I know cos(alpha) from the dot product // sin(theta) = -MatrixDot(Ghat,kihat) Variable lambda = 4*PI/G * (-MatrixDot(Ghat,kihat)) // 4 sin(theta) / G lambda = (lambda<0) ? NaN : lambda if (printIt) printf " [%s] x (%d %d %d) --> lambda=%g\r",NameOfWave(recip),h,k,l,lambda endif KillWaves/Z Ghat, kihat, hkl__ return lambda End // find where (hkl) should land on the detector from reciprocal lattice, and print to history Function/T printDetectorXY_of_hkl(recip,h,k,l) Wave/D recip // 3x3 matrix with reciprocal lattice Variable h,k,l if (numtype(h) || numtype(k) || numtype(l) || !WaveExists(recip)) h = numtype(h) ? 0 : h k = numtype(k) ? 2 : k l = numtype(l) ? 2 : l Prompt h, "h" Prompt k, "k" Prompt l, "l" String recipLattice, mat3x3List=WaveList_NxN(3,3) recipLattice = SelectString(exists("recip")==1,"","recip") Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List DoPrompt "hkl",recipLattice,h,k,l if (V_flag) return "" endif Wave/D recip = $recipLattice endif if (numtype(h) || numtype(k) || numtype(l) || !WaveExists(recip)) return "" endif String graphName = StringFromList(0,WinList("*",";","WIN:")) if (strlen(graphName)<1) Abort "No Window for doing determining parameters" endif Wave image = ImageNameToWaveRef(graphName,StringFromList(0,ImageNameList(graphName,";"))) if (!WaveExists(image)) return "" endif Variable detX, detY // position on detector (mm) Make/O hkl__temp__ = {h,k,l} Wave hkl = hkl__temp__ calcXYofHKL(recip,hkl,detX,detY,note(image)) // find xy on the detector printf "for hkl=(%d %d %d) [%s] --> detector=(%g, %g)\r",hkl[0],hkl[1],hkl[2],NameOfWave(recip),detX,detY KillWaves hkl__temp__ String result sprintf result "detX=%g;detY=%g;",detX,detY return result End Function/T hkl_BeamAxis(recip) // computes hkl of a Gvector that is from sample back toward the source Wave/D recip Variable printIt=0 if (!WaveExists(recip) || DimSize(recip,0)!=3 || DimSize(recip,1)!=3 || DimSize(recip,2)!=0) String recipLattice, mat3x3List=WaveList_NxN(3,3) recipLattice = SelectString(exists("recip")==1,"","recip") Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List DoPrompt "reciprocal lattice",recipLattice if (V_flag) return "" endif Wave/D recip = $recipLattice printIt = 1 endif // recip x H = G --> recipInv x recip x H = recipInv x G // hkl = recipInv x Gvec, so find inverse of recip Make/N=3/O Gvec={0,0,-2*PI} // back reflection for lambda=1 Make/N=(3,3)/O/D recip_temp_,recipInv recipInv = (p==q) recip_temp_ = recip MatrixLinearSolve /O recip_temp_ recipInv // recipInv now contains the inverse of hklN KillWaves/Z recip_temp_ MatrixMultiply recipInv, Gvec // result is in M_product Wave hkl = M_product String out sprintf out, "%.2f, %.2f, %.2f",hkl[0],hkl[1],hkl[2] KillWaves/Z Gvec, W_IPIV, recipInv, M_product if (printIt) printf "[%s]x[hkl=(%s)] = G, points at the source\r",NameOfWave(recip),out endif return out End Function addPredictedHKL(recip, hklMax,lamMin) Wave/D recip Variable hklMax Variable lamMin // usually 0.42 if (!WaveExists(recip) || numtype(hklMax) || hklMax<1 || lamMin<0 || numtype(lamMin)) lamMin = (lamMin<0 || numtype(lamMin)) ? 0.42 : lamMin String mat3x3List=WaveList_NxN(3,3) String recipLattice = SelectString(exists("recip")==1,"","recip") hklMax = (numtype(hklMax) || hklMax<1)? 6 : hklMax Prompt hklMax, "max value for h,k, or l" Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List Prompt lamMin, "minimum allowed lambda ()" DoPrompt "choose max", recipLattice, hklMax,lamMin Wave/D recip = $recipLattice if (V_flag || !WaveExists(recip)) return NaN endif printf " addPredictedHKL(%s,%d,%g)\r",recipLattice,hklMax,lamMin endif SVAR materialsList=root:Packages:SCD:materialsList String material = StringByKey("material",note(recip),"=") Variable i=WhichListItem(material,materialsList) if (i<0) Abort "unknown material '"+material+"'" endif SVAR BravaisList=root:Packages:SCD:BravaisList String funcName = "isLowestOrder"+StringFromList(i,BravaisList) FUNCREF isLowestOrderFCC isLowestOrder = $funcName // function to use to determine lowest hkl Variable h,k,l Make/N=3/O/D hkl SetDrawLayer /K UserFront // clear the drawing layer for (l=0;l<=hklMax;l+=(l>=0)) l = -l for (k=0;k<=hklMax;k+=(k>=0)) k = -k for (h=0;h<=hklMax;h+=(h>=0)) h = -h hkl = {h,k,l} if (isLowestOrder(hkl)) placeOneSpot(recip,hkl,lamMin) endif endfor endfor endfor KillWaves/Z hkl End //Function addPredictedHKL(recip, hklMax,lamMin) // Wave/D recip // Variable hklMax // Variable lamMin // usually 0.42 // if (!WaveExists(recip) || numtype(hklMax) || hklMax<1 || lamMin<0 || numtype(lamMin)) // lamMin = (lamMin<0 || numtype(lamMin)) ? 0.42 : lamMin // String mat3x3List=WaveList_NxN(3,3) // String recipLattice = SelectString(exists("recip")==1,"","recip") // hklMax = (numtype(hklMax) || hklMax<1)? 6 : hklMax // Prompt hklMax, "max value for h,k, or l" // Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List // Prompt lamMin, "minimum allowed lambda ()" // DoPrompt "choose max", recipLattice, hklMax,lamMin // Wave/D recip = $recipLattice // if (V_flag || !WaveExists(recip)) // return NaN // endif // printf " addPredictedHKL(%s,%d,%g)\r",recipLattice,hklMax,lamMin // endif // // SVAR materialsList=root:Packages:SCD:materialsList // String material = StringByKey("material",note(recip),"=") // Variable i=WhichListItem(material,materialsList) // if (i<0) // Abort "unknown material '"+material+"'" // endif // SVAR BravaisList=root:Packages:SCD:BravaisList // String funcName = "isLowestOrder"+StringFromList(i,BravaisList) // FUNCREF isLowestOrderFCC isLowestOrder = $funcName // function to use to determine lowest hkl // // Variable h,k,l // Make/N=3/O/D hkl // SetDrawLayer /K UserFront // clear the drawing layer // // for (l=-hklMax;l<=hklMax;l+=1) // for (k=-hklMax;k<=hklMax;k+=1) // for (h=-hklMax;h<=hklMax;h+=1) // hkl = {h,k,l} // if (isLowestOrder(hkl)) // placeOneSpot(recip,hkl,lamMin) // endif // endfor // endfor // endfor // KillWaves/Z hkl //End Function placeOneSpot(recip,hkl,lamMin) // place hkl spots on top image according to reciprocal lattice vectors // coordinate system is the standard beamline coordinate system, z-along beam, y-up, x-sideways Wave/D recip // a 3x3 reciprocal lattice a*=recip[][0], b*=recip[][1], c*=recip[][2] Wave/D hkl // hkl vector Variable lamMin // minimum allowed wavelength Variable lambda = hkl2lambda(recip,hkl[0],hkl[1],hkl[2]) if (lambdadetX || detX>xhi || ylo>detY || detY>yhi) // spot not on detector return NaN endif printf "drawing for (%2d %2d %2d) at %4.2f () at detector=(%.2f,%.2f)\r",hkl[0],hkl[1],hkl[2],lambda,detX,detY Variable boxSize=3 SetDrawLayer UserFront // SetDrawEnv xcoord= bottom,ycoord= left,fillpat= 0 SetDrawEnv xcoord= bottom,ycoord= left,linepat= 3,fillpat= 0 DrawRect detX-boxSize,detY+boxSize, detX+boxSize,detY-boxSize // DrawRect left, top, right, bottom String str if (detY>-100) SetDrawEnv xcoord=bottom,ycoord=left,fsize=9,textxjust=1,textyjust=2 else SetDrawEnv xcoord=bottom,ycoord=left,fsize=9,textxjust=1,textyjust=0 endif sprintf str, "(%d %d %d)",hkl[0],hkl[1],hkl[2] DrawText detX,detY,str // DrawText x0, y0, textStr End Function calcXYofHKL(recip,hkl,detX,detY,noteStr) // calc place of hkl on detector from reciprocal lattice vectors // coordinate system is the standard beamline coordinate system, z-along beam, y-up, x-sideways Wave/D recip // a 3x3 reciprocal lattice a*=recip[][0], b*=recip[][1], c*=recip[][2] Wave/D hkl // hkl vector Variable &detX, &detY // position on detector (mm), passed back String noteStr detX = NaN ; detY = NaN // default is an error Variable G // length of G vector MatrixMultiply recip, hkl // result is in M_product KillWaves/Z Gvec Rename M_product Gvec // G vector (from hkl) Wave Gvec=Gvec Redimension/N=3 Gvec Duplicate/O/D Gvec Ghat G = normalize(Ghat) // Ghat is now unit vector in Gvec direction, and G = |Gvec| Make/O/D/N=3 kihat={0,0,1}, kfhat, detCentHat, inDetPlane, detXhat Variable k = -G/2 / MatrixDot(Ghat,kihat) // = 2/lambda, -ki^ G^ = sin(theta) // to solve this remember that kf = ki + Gvec kfhat = (k*kihat[p]) + Gvec[p] // actually I wil only need kfhat, not kf normalize(kfhat) Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta") Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd") Variable cosDeta=cos(deta*PI/180), sinDeta=sin(deta*PI/180) detCentHat = {sinDeta, 0, cosDeta} // unit vector from (000) to center of detector Variable dot = MatrixDot(detCentHat,kfhat) if (dot<0.01) return NaN // kfhat is does not go towards detector (parallel or worse) endif Variable d = detd /dot // distance from origin to detector of scattered ray (mm) // vector from detector center to pixel = (origin to pixel) - (origin to detector center) inDetPlane = d*kfhat[p] - detd*detCentHat[p] // detd/dot is distance from origin pixel on detector (mm) detXhat = {-cosDeta,0,sinDeta} // unit vector from detector center toward +x on detector detX = MatrixDot(inDetPlane,detXhat) detY = inDetPlane[1] KillWaves/Z Gvec, Ghat, kihat, kfhat, detCentHat, inDetPlane, detXhat return d End Function spots2recip(x0,y0,x1,y1) // convert 2 spots on detector to reciprocal lattice Variable x0,y0,x1,y1 x0 = numtype(x0) ? hcsr(A) : x0 // position on detector (mm), in detector frame y0 = numtype(y0) ? vcsr(A) : y0 x1 = numtype(x1) ? hcsr(B) : x1 y1 = numtype(y1) ? vcsr(B) : y1 Make/N=3/O/D hkl0,hkl1, hkl2 if (numtype(x0) || numtype(y0) || numtype(x1) || numtype(y1)) Abort "need cursors on spots" endif SVAR materialsList=root:Packages:SCD:materialsList String material Prompt material, "material to index", popup, materialsList DoPrompt "pick material", material if (V_Flag) return 1 endif Variable i=WhichListItem(material,materialsList) if (i<0) Abort "unknown material '"+material+"'" endif SVAR aoList=root:Packages:SCD:aoList Variable ao = str2num(StringFromList(i,aoList)) // d-spaceing () Variable h,k,l Prompt h, "h" Prompt k, "k" Prompt l, "l" h=0 ; k=2 ; l=2 DoPrompt "hkl of first point (cursor A)",h,k,l if (V_flag) return NaN endif hkl0 = {h,k,l} h=1 ; k=3 ; l=3 DoPrompt "hkl of second point (cursor B)",h,k,l if (V_flag) return NaN endif hkl1 = {h,k,l} String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";")))) String recipName sprintf recipName, "recip_%s_%d", material, 10*NumberByKey("depth",noteStr,"=") // I have all of the info I need here, construct the reciprocal lattice. Make/O/D/N=3 kihat={0,0,1}, kfhat, G0,G1,G2 // kihat points down stream, from source to sample Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta") Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd") Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx") Variable xp, zp // pixel position (in horizontal plane) pixelXZfromDetector(deta,detd,x0,xp,zp) // calculate pixel coordinates xp,zp // kfhat = {xp,y0,zp} // vertical height (y-direction) is y0 kfhat = {xp-bx,y0,zp} // vertical height (y-direction) is y0 normalize(kfhat) G0 = kfhat-kihat // points half way between kf and -ki, normalize later pixelXZfromDetector(deta,detd,x1,xp,zp) // calculate pixel coordinates xp,zp // kfhat = {xp,y1,zp} // vertical height (y-direction) is y1 kfhat = {xp-bx,y1,zp} // vertical height (y-direction) is y1 normalize(kfhat) G1 = kfhat-kihat crossV4(G0,G1,G2) // vector cross product to get third direction crossV4(hkl0,hkl1,hkl2) RemoveCommonFactors(hkl0) RemoveCommonFactors(hkl1) RemoveCommonFactors(hkl2) printwave(hkl0) printwave(hkl1) printwave(hkl2) printwave(G0) printwave(G1) printwave(G2) printf "creating the reciprocal lattice in wave '%s'\r",recipName print "" crossV4(G2,G0,G1) // ensure an orthonormal system crossV4(hkl2,hkl0,hkl1) doAlert 1, "do the funky 'G2 *= -1'?" if(V_Flag==1) G2 *= -1 // this seems to be necessary, but I do not know why !!!!!! print " did the funky 'G2 *= -1'" else print " took out the 'G2 *= -1'" endif normalize(G0) normalize(G1) normalize(G2) normalize(hkl0) normalize(hkl1) normalize(hkl2) printwave(hkl0) printwave(hkl1) printwave(hkl2) printwave(G0) printwave(G1) printwave(G2) // now G0,G1,G2 and hkl0,hkl1,hkl2 are all defined as orthonormal systems // form a matrix of the hklN // recip x H = G --> recip = G x H^-1 Make/N=(3,3)/D/O hklN, ident, $recipName Wave/D recip = $recipName noteStr = ReplaceStringByKey("material",noteStr,material,"=") Note/K recip Note recip noteStr ident = (p==q) // identity matrix hklN[][0] = hkl0[p] hklN[][1] = hkl1[p] hklN[][2] = hkl2[p] print "determinant(hklN)= ",MatrixDet(hklN) MatrixLinearSolve /O hklN ident // ident now contains the inverse of hklN KillWaves/Z W_IPIV recip[][0] = G0[p] recip[][1] = G1[p] recip[][2] = G2[p] print "determinant(G012)= ",MatrixDet(recip) MatrixMultiply recip, ident // multiplying G x H^-1 Wave M_product=M_product recip = M_product * (2*PI/ao) print "determinant(recip)= ",MatrixDet(recip) hkl0 = {0, 1, 2} hkl1 = {1, 2, 2} hkl2 = {-2, 2, -1} print "" MatrixMultiply recip, hkl0 printWave(hkl0) normalize(M_product) printWave(M_product) MatrixMultiply recip, hkl1 printWave(hkl1) normalize(M_product) M_product *= 1.60552 printWave(M_product) MatrixMultiply recip, hkl2 printWave(hkl2) normalize(M_product) printWave(M_product) KillWaves/Z G0,G1,G2, hkl0,hkl1,hkl2, M_product,kfhat,kihat, hklN, ident End // supplanted by the Cross operation in Igor 5 Function crossV4(a,b,c) // vector cross product Wave a,b,c c[2] = a[0]*b[1] - a[1]*b[0] c[0] = a[1]*b[2] - a[2]*b[1] c[1] = a[2]*b[0] - a[0]*b[2] End // ==================================================================== // ========================== End of Reciprocal Lattice ======================= // ==================================================================== // ==================================================================== // ==================================================================== // ==================================================================== // ============================= Start of Movies ========================== // ==================================================================== Function RunAsMovie() GetWindow $StringFromList(0,WinList("*",";","WIN:")), note Variable oldVal = NumberByKey("value",S_value,"=") upDateSliceProc("",-inf,1) // goto the start DoUpdate Sleep/T 30 Variable i = 0 do Sleep /T 4 i = upDateSliceProc("up",1,1) ? i+1 : i DoUpdate while (i<2) Sleep /T 30 upDateSliceProc("",oldVal,1) // goto the original position End Function makeMovieOfTopSliceGraph(upDown) //each frame is at one step in the slice graph Variable upDown // direction of movie, true = up, false=down if (numtype(upDown)) upDown = 1 Prompt upDown, "direction", popup,"up;down" DoPrompt "direction of movie", upDown if (V_flag) return 1 endif upDown = (upDown==1) endif String winNam=StringFromList(0,WinList("*",";","WIN:")) if (strlen(winNam)<1) DoAlert 0, "No Window to make movie of" return 1 endif GetWindow $winNam, note String noteStr = S_value if (NumberByKey("sliceWindow",noteStr,"=")!=1) DoAlert 0, "There is no slice window as top graph" return 1 endif if (upDown) upDateSliceProc("",-inf,1) // goto the start else upDateSliceProc("",inf,1) // goto the start endif String source = StringByKey("source",noteStr,"=") printf "Creating a movie named '%s.mov'\r",source NewMovie /F=10/P=home/Z as source+".mov" if (V_flag) DoAlert 0, "Could not create movie" return 1 endif Variable i = 0 do DoUpdate AddMovieFrame if (upDown) i = upDateSliceProc("up",1,1) ? i+1 : i else i = upDateSliceProc("down",1,1) ? i+1 : i endif while (i<2) CloseMovie End // ==================================================================== // ============================== End of Movies ========================== // ==================================================================== // ==================================================================== // ==================================================================== // ==================================================================== // ============================= Start of hkl annotations ===================== // ==================================================================== // this form is to be called from a marquee around a spot Function EvaluateSpot(): GraphMarquee GetMarquee/K left, bottom if (!V_flag) return NaN endif Wave image = ImageNameToWaveRef(S_marqueeWin,StringFromList(0,ImageNameList(S_marqueeWin,";"))) if (!WaveExists(image)) return NaN endif // String format = "flag: %g; left: %g; top: %g; right: %g; bottom: %g window:'%s'\r" // printf format, V_flag, V_left, V_top, V_right, V_bottom, S_marqueeWin Make/O/T/N=2 T_Constraints T_Constraints[0] = {"K6 > 0","K6 < 1"} Make/D/O W_coef={0,0,0,0,0,0,0.5} Variable V_FitError, V_fitOptions=4 CurveFit/Q/N Gauss2D kwCWave=W_coef, image(V_left,V_right)(V_bottom,V_top) /C=T_Constraints KillWaves/Z T_Constraints if (V_FitError) DoAlert 0, "could not fit peak" return NaN endif Variable tolerance=.13 // tolerated error in sqrt(h^2 + k^2 + l^2) Variable detX = w_coef[2] // position on detector (mm) Variable detY = w_coef[4] Variable tth = calc2theta(detX,detY) Variable lambda = getLambdaFromPlot("") Variable Integral = integralGauss2d(V_left,V_right,V_bottom,V_top,w_coef[0],image) if (numtype(lambda) || tth>=PI || tth<0) // no wavelength, just get 2theta // printf "at detector=(%.3f,%.3f), 2theta = %.4fr",detX,detY,tth*180/PI printf "at detector=(%.3f,%.3f), 2theta = %.4f, Integral = %g\r",detX,detY,tth*180/PI,Integral return 1 endif KillWaves/Z W_coef,W_sigma,W_ParamConfidenceInterval SVAR materialsList=root:Packages:SCD:materialsList String material="all" if (ItemsInList(materialsList)>1) Prompt material, "material to index", popup, materialsList+";-;all" DoPrompt "material", material if (V_flag) return 1 endif endif String str,result Variable i if (!cmpstr("all",material)) // check all the materials in the materialsList Variable m,N = ItemsInList(materialsList) String mat for (m=0;m0) text += "\r" endif text += StringByKey("hkl",item,"=",",") endif endfor DoAlert 1, "add label to plot?" if (V_flag==1) i = WhichListItem(material,materialsList) str = "\\K("+SelectString(i-1,"65535,0,0","0,0,65535","0,65535,0")+")" addAnnotaionHKL(S_marqueeWin,detX,detY,str+text) endif if (exists("recip")) Wave/D recip=recip printf "using reciprocal lattice, hkl = (%s)\r",hkl_fromDetectorSpot_recip(detX,detY,recip) endif End Function integralGauss2d(V_left,V_right,V_bottom,V_top,K0,image) Variable V_left,V_right,V_bottom,V_top // size of box in axis units (mm) Variable K0 Wave image Variable ilo = (V_left-DimOffset(image,0)) / DimDelta(image,0) Variable ihi = (V_right-DimOffset(image,0)) / DimDelta(image,0) Variable jlo = (V_bottom-DimOffset(image,1)) / DimDelta(image,1) Variable jhi = (V_top-DimOffset(image,1)) / DimDelta(image,1) Variable i,j, N=0, Int=0 for (j=jlo;j<=jhi;j+=1) for (i=ilo;i<=ihi;i+=1) Int += image[i][j] N += 1 endfor endfor Int -= N*K0 // remove background Int *= DimDelta(image,0)* DimDelta(image,1) // scale by dx,dy to make it an area return Int // K0+K1*exp((-1/(2*(1-K6^2)))*(((x-K2)/K3)^2 + ((y-K4)/K5)^2 - (2*K6*(x-K2)*(y-K4)/(K3*K5)))) End // this form is to be called from a menu, or command line Function InfoAboutSpotOnDetector(detX,detY,tolerance) Variable detX, detY // postion on detector (mm), or NaN to use cursor A Variable tolerance // tolerated error in sqrt(h^2 + k^2 + l^2) (=.1) Variable tth = calc2theta(detX,detY) printf "at detector=(%.3f,%.3f), 2theta = %.4fr",detX,detY,tth*180/PI if (numtype(getLambdaFromPlot(""))) // good wavelength, try for an hkl return tth endif SVAR materialsList=root:Packages:SCD:materialsList String material Prompt material, "material to index", popup, materialsList DoPrompt "material", material if (V_flag) return 1 endif String str = calcHKL(detX,detY,material,tolerance) if (strlen(str)<1) print "no hkl close enough to point on detector" return tth endif Variable i for (i=0;i0 || DimSize(v2,1)>0 || DimSize(v1,0)!=DimSize(v2,0)) // 1-d vectors must match // return NaN // endif // Variable N = DimSize(v1,0) // Variable i,m1=0,m2=0, dot=0 // for (i=0;i 1 // return acos(dot) //End // =============================================================== Function addAnnotaionHKL(graphName,xx,yy,str) String graphName Variable xx,yy // postion of annotaion String str // string to add String list = AnnotationList(graphName) String name // annotaion name Variable i=0 do name = "hkl"+num2istr(i) i += 1 while(WhichListItem(name,list)>0) // this loop finds next name of type hklN Wave image = ImageNameToWaveRef(graphName,StringFromList(0,ImageNameList(graphName,";"))) if (!WaveExists(image)) return NaN endif Variable index // index into image corresponding to xx,yy index = round((yy-DimOffset(image,1))/DimDelta(image,1)) // index in y direction index *= DimSize(image,0) // ponits to correct row index += round((xx-DimOffset(image,0))/DimDelta(image,0)) // along row index = pnt2x(image,index) // funny thing Igor needs // str = "\\Z09\\K(65535,0,0)"+str str = "\\Z09"+str Tag/W=$graphName/A=MC/B=1/F=0/X=0/Y=0/N=$name/I=1/L=0/Z=1 $NameOfWave(image),index,str End Function deleteAnnotationsHKL(): GraphMarquee GetMarquee/K left, bottom if (!V_flag) return NaN endif String list = AnnotationList(S_marqueeWin) String name Variable i, N=ItemsInList(list) for (i=0;i= thresh, 0 includes all // Variable res=1.036/2 // resolution half-width (mm), assumed flat topped Variable printIt=0 if (numtype(depth+tol+thresh) || tol<=0 || thresh<0 || strlen(base)<1 || strlen(range)<1) depth = numtype(depth) ? 0 : depth tol = (numtype(tol) || tol<=0) ? 0.5 : tol thresh = numtype(thresh) ? 0 : abs(thresh) base = SelectString(strlen(base) , "stackA", base) SVAR rangeList = root:Packages:SCD:rangeList range = SelectString(strlen(range),StringFromList(ItemsInList(rangeList)-1,rangeList),range) Prompt depth, "desired depth in the sample (mm)" Prompt tol, "tolerance on the depth (depth眛ol) (mm)" Prompt thresh, "threshold for including a pixel, 0 includes all" Prompt base, "for wave name, e.g. base='stackA'" Prompt range, "range of run numbers in original data set, e.g. '8523-8538'",popup,rangeList+";other" DoPrompt "the depth to sort out", depth, tol,base,range,thresh if (V_flag) return "" endif if (!cmpstr(range,"other")) Prompt range, "range of run numbers, e.g. '3-10,12'" DoPrompt "new range", range if (V_flag) return "" endif rangeList = AddListItem(range,rangeList,";",Inf) endif printIt = 1 endif thresh = round(thresh) if (numtype(depth+tol+thresh) || thresh<0 || tol<=0 || strlen(base)<1 || strlen(range)<1) print "Illegal depth or tolerance" return "" endif String sliceName sprintf sliceName "%s_%s%02d_%02d",base, SelectString(depth<0,"","m"),abs(depth*10),10*tol if (printIt) printf "SortDataOneDepth(%g,%g,\"%s\",\"%s\",%g)\t\t\t// into: '%s'\r",depth,tol,base,range,thresh,sliceName endif range = expandRange(range,";") String frameName sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(0,range)) Wave frame = $frameName Variable dimX=DimSize(frame,0),dimY=DimSize(frame,1),dimZ=DimSize(frame,2) String noteStr = note(frame) Make/N=(dimX,dimY,dimZ)/O/U/W $sliceName // create 1 depth slice for x, y, and lambda Wave slice = $sliceName slice = 0 SetXYtoANGER(slice) noteStr = ReplaceNumberByKey("depth",noteStr,depth,"=") noteStr = ReplaceNumberByKey("tol",noteStr,tol,"=") noteStr = ReplaceNumberByKey("threshSort",noteStr,thresh,"=") noteStr = RemoveByKey("translate_mm",noteStr,"=") Note slice, noteStr Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta") Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd") Variable x2mm = valFromWaveNoteOrGlobal(noteStr,"x2mm") Variable y2mm = valFromWaveNoteOrGlobal(noteStr,"y2mm") Variable xLeftmm = valFromWaveNoteOrGlobal(noteStr,"xLeftmm") Variable slitx = valFromWaveNoteOrGlobal(noteStr,"slitx") Variable slitz = valFromWaveNoteOrGlobal(noteStr,"slitz") Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx") Variable surface = valFromWaveNoteOrGlobal(noteStr,"surface") Variable k // index over each of the runs in range Variable i // index over each pixels in the x-direcion in the detector Variable trans // translation position for each SCD run Variable xi // x position on detector from its center (mm) Variable xp, zp // x,z coodinates of the pixel (mm), in beam line system Variable zintercept // position on beam axis where beam hits (mm) // Variable x0=DimOffset(slice,0) // Variable dx=DimDelta(slice,0) // Variable dmin, dmax // range of the overlap Variable/G V_nonZero=0 // flags if this depth is empty Variable angle // angle between detector normal and exit ray Variable detaR = deta*PI/180 // radian version of deta Variable iorigin // index where beam would land if this volume were at (0,0) Variable start = ticks do // look at each run in the range if (strlen(StringFromList(k,range))<1) break endif sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(k,range)) Wave frame = $frameName if (!WaveExists(frame)) Abort "the input wave '"+frameName+"' cannot be found" endif trans = NumberByKey("translate_mm",note(frame),"=") for (i=0;iPI/2) printf "angle = %g-%g = %g, i = %d, iorigin = %d, %s\r",deta,atan2((xp-bx),(zp-zintercept))*180/PI,angle*180/PI,i,iorigin,NameOfWave(frame) elseif (thresh<0.5) slice[iorigin][][] += frame[i][q][r] else slice[iorigin][][] += (frame[i][q][r]>thresh) ? frame[i][q][r] : 0 endif endif // dmin = max(zmm-depth-res,depth-tol) // bottom of overlap range // dmax = min(zmm-depth+res,depth+tol) // top of overlap range // if (dmax>dmin) // slice[i][][] += round((dmax-dmin)/res/2*frame[i][q][r]) // endif endfor k += 1 while(1) if (printIt) printf "sorting took %.2f sec\r",(ticks-start)/60.15 DoAlert 1, "put up a plot of "+sliceName if (V_flag) ShowSliceGraph(sliceName) endif endif return GetWavesDataFolder(slice,2) // return the full path End Function/T SortDataOneDepthOLD(depth,tol,base,range,thresh) // sort data into 1 depth array // sort data and save only one depth at the given depth (the result is 3-d, x y & lambda) // returns the full path name of the new 3-d wave // Note that +z is along beam, +y is up, and +x is horizontally perp to beam (points to where you stand). Variable depth // target depth to get (mm) Variable tol // tolerance on the depth (mm), typically 0.5 String base // for file name, e.g. base="CuSCD", or "stackA" String range // range of run numbers in original data set Variable thresh // only include pixels that are >= thresh, 0 includes all // Variable res=1.036/2 // resolution half-width (mm), assumed flat topped Variable printIt=0 if (numtype(depth+tol+thresh) || tol<=0 || thresh<0 || strlen(base)<1 || strlen(range)<1) depth = numtype(depth) ? 0 : depth tol = (numtype(tol) || tol<=0) ? 0.5 : tol thresh = numtype(thresh) ? 0 : abs(thresh) base = SelectString(strlen(base) , "stackA", base) SVAR rangeList = root:Packages:SCD:rangeList range = SelectString(strlen(range),StringFromList(ItemsInList(rangeList)-1,rangeList),range) Prompt depth, "desired depth in the sample (mm)" Prompt tol, "tolerance on the depth (depth眛ol) (mm)" Prompt thresh, "threshold for including a pixel, 0 includes all" Prompt base, "for wave name, e.g. base='stackA'" Prompt range, "range of run numbers in original data set, e.g. '8523-8538'",popup,rangeList+";other" DoPrompt "the depth to sort out", depth, tol,base,range,thresh if (V_flag) return "" endif if (!cmpstr(range,"other")) Prompt range, "range of run numbers, e.g. '3-10,12'" DoPrompt "new range", range if (V_flag) return "" endif rangeList = AddListItem(range,rangeList,";",Inf) endif printIt = 1 endif thresh = round(thresh) if (numtype(depth+tol+thresh) || thresh<0 || tol<=0 || strlen(base)<1 || strlen(range)<1) print "Illegal depth or tolerance" return "" endif String sliceName sprintf sliceName "%s_%s%02d_%02d",base, SelectString(depth<0,"","m"),abs(depth*10),10*tol if (printIt) printf "SortDataOneDepth(%g,%g,\"%s\",\"%s\",%g)\t\t\t// into: '%s'\r",depth,tol,base,range,thresh,sliceName endif range = expandRange(range,";") String frameName sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(0,range)) Wave frame = $frameName Variable dimX=DimSize(frame,0),dimY=DimSize(frame,1),dimZ=DimSize(frame,2) String noteStr = note(frame) Make/N=(dimX,dimY,dimZ)/O/U/W $sliceName // create 1 depth slice for x, y, and lambda Wave slice = $sliceName slice = 0 SetXYtoANGER(slice) noteStr = ReplaceNumberByKey("depth",noteStr,depth,"=") noteStr = ReplaceNumberByKey("tol",noteStr,tol,"=") noteStr = ReplaceNumberByKey("threshSort",noteStr,thresh,"=") noteStr = RemoveByKey("translate_mm",noteStr,"=") Note slice, noteStr Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta") Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd") Variable x2mm = valFromWaveNoteOrGlobal(noteStr,"x2mm") Variable y2mm = valFromWaveNoteOrGlobal(noteStr,"y2mm") Variable xLeftmm = valFromWaveNoteOrGlobal(noteStr,"xLeftmm") Variable slitx = valFromWaveNoteOrGlobal(noteStr,"slitx") Variable slitz = valFromWaveNoteOrGlobal(noteStr,"slitz") Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx") Variable surface = valFromWaveNoteOrGlobal(noteStr,"surface") Variable k // index over each of the runs in range Variable i // index over each pixels in the x-direcion in the detector Variable trans // translation position for each SCD run Variable xi // x position on detector from its center (mm) Variable xp, zp // x,z coodinates of the pixel (mm), in beam line system Variable zintercept // position on beam axis where beam hits (mm) // Variable x0=DimOffset(slice,0) // Variable dx=DimDelta(slice,0) // Variable dmin, dmax // range of the overlap Variable/G V_nonZero=0 // flags if this depth is empty Variable start = ticks do // look at each run in the range if (strlen(StringFromList(k,range))<1) break endif sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(k,range)) Wave frame = $frameName if (!WaveExists(frame)) Abort "the input wave '"+frameName+"' cannot be found" endif trans = NumberByKey("translate_mm",note(frame),"=") for (i=0;ithresh) ? frame[i][q][r] : 0 endif endif // dmin = max(zmm-depth-res,depth-tol) // bottom of overlap range // dmax = min(zmm-depth+res,depth+tol) // top of overlap range // if (dmax>dmin) // slice[i][][] += round((dmax-dmin)/res/2*frame[i][q][r]) // endif endfor k += 1 while(1) if (printIt) printf "sorting took %.2f sec\r",(ticks-start)/60.15 DoAlert 1, "put up a plot of "+sliceName if (V_flag) ShowSliceGraph(sliceName) endif endif return GetWavesDataFolder(slice,2) // return the full path End Function/T SortDataOneDepth_RealOLD(depth,tol,base,range,thresh) // sort data into 1 depth array // sort data and save only one depth at the given depth (the result is 3-d, x y & lambda) // returns the full path name of the new 3-d wave // Note that +z is along beam, +y is up, and +x is horizontally perp to beam (points to where you stand). Variable depth // target depth to get (mm) Variable tol // tolerance on the depth (mm), typically 0.5 String base // for file name, e.g. base="CuSCD", or "stackA" String range // range of run numbers in original data set Variable thresh // only include pixels that are >= thresh, 0 includes all Variable res=1.036/2 // resolution half-width (mm), assumed flat topped Variable printIt=0 if (numtype(depth+tol+thresh) || tol<=0 || thresh<0 || strlen(base)<1 || strlen(range)<1) depth = numtype(depth) ? 0 : depth tol = (numtype(tol) || tol<=0) ? 0.5 : tol thresh = numtype(thresh) ? 0 : abs(thresh) base = SelectString(strlen(base) , "stackA", base) SVAR rangeList = root:Packages:SCD:rangeList range = SelectString(strlen(range),StringFromList(ItemsInList(rangeList)-1,rangeList),range) Prompt depth, "desired depth in the sample (mm)" Prompt tol, "tolerance on the depth (depth眛ol) (mm)" Prompt thresh, "threshold for including a pixel, 0 includes all" Prompt base, "for wave name, e.g. base='stackA'" Prompt range, "range of run numbers in original data set, e.g. '8523-8538'",popup,rangeList+";other" DoPrompt "the depth to sort out", depth, tol,base,range,thresh if (V_flag) return "" endif if (!cmpstr(range,"other")) Prompt range, "range of run numbers, e.g. '3-10,12'" DoPrompt "new range", range if (V_flag) return "" endif rangeList = AddListItem(range,rangeList,";",Inf) endif printIt = 1 endif thresh = round(thresh) if (numtype(depth+tol+thresh) || thresh<0 || tol<=0 || strlen(base)<1 || strlen(range)<1) print "Illegal depth or tolerance" return "" endif String sliceName sprintf sliceName "%s_%s%02d_%02d",base, SelectString(depth<0,"","m"),abs(depth*10),10*tol if (printIt) printf "SortDataOneDepth(%g,%g,\"%s\",\"%s\",%g)\t\t\t// into: '%s'\r",depth,tol,base,range,thresh,sliceName endif range = expandRange(range,";") String frameName sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(0,range)) Wave frame = $frameName Variable dimX=DimSize(frame,0),dimY=DimSize(frame,1),dimZ=DimSize(frame,2) String noteStr = note(frame) Make/N=(dimX,dimY,dimZ)/O/U/W $sliceName // create 1 depth slice for x, y, and lambda Wave slice = $sliceName slice = 0 SetXYtoANGER(slice) noteStr = ReplaceNumberByKey("depth",noteStr,depth,"=") noteStr = ReplaceNumberByKey("tol",noteStr,tol,"=") noteStr = ReplaceNumberByKey("threshSort",noteStr,tol,"=") Note slice, noteStr Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta") Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd") Variable x2mm = valFromWaveNoteOrGlobal(noteStr,"x2mm") Variable y2mm = valFromWaveNoteOrGlobal(noteStr,"y2mm") Variable xLeftmm = valFromWaveNoteOrGlobal(noteStr,"xLeftmm") Variable slitx = valFromWaveNoteOrGlobal(noteStr,"slitx") Variable slitz = valFromWaveNoteOrGlobal(noteStr,"slitz") Variable surface = valFromWaveNoteOrGlobal(noteStr,"surface") Variable k // index over each of the runs in range Variable i // index over each pixels in the x-direcion in the detector Variable trans // translation position for each SCD run Variable xi // x position on detector from its center (mm) Variable xp, zp // x,z coodinates of the pixel (mm), in beam line system Variable zintercept // position on beam axis where beam hits (mm) // Variable x0=DimOffset(slice,0) // Variable dx=DimDelta(slice,0) // Variable dmin, dmax // range of the overlap Variable start = ticks do // look at each run in the range if (strlen(StringFromList(k,range))<1) break endif sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(k,range)) Wave frame = $frameName if (!WaveExists(frame)) Abort "the input wave '"+frameName+"' does cannot be found" endif trans = NumberByKey("translate_mm",note(frame),"=") for (i=0;ithresh) ? frame[i][q][r] : 0 endif endif // dmin = max(zmm-depth-res,depth-tol) // bottom of overlap range // dmax = min(zmm-depth+res,depth+tol) // top of overlap range // if (dmax>dmin) // slice[i][][] += round((dmax-dmin)/res/2*frame[i][q][r]) // endif endfor k += 1 while(1) printf "sorting took %.2f sec\r",(ticks-start)/60.15 if (printIt) DoAlert 1, "put up a plot of "+sliceName if (V_flag) ShowSliceGraph(sliceName) endif endif return GetWavesDataFolder(slice,2) // return the full path End // this routine extracts the value of varName from noteStr if it is there, otherwise, it tries to // use the global the global of same name in root:Packages:SCD. If all fail, it returns NaN. Function valFromWaveNoteOrGlobal(noteStr,varName) String noteStr String varName Variable val = NumberByKey(varName,noteStr,"=") if (numtype(val)) String gName = "root:Packages:SCD:"+varName val = NumVarOrDefault(gName,NaN) endif return val End // ==================================================================== // =========================== End Sorting and Summing ===================== // ==================================================================== // ==================================================================== // ==================================================================== // ==================================================================== // ============================== Start Displaying ========================= // ==================================================================== Function ShowSliceGraph(source) String source if (exists(source)!=1) source = "" elseif (DimSize($source,2)<1) source = "" endif if (strlen(source)<1) source = StringFromList(ItemsInList(WaveList("*",";","DIMS:3"))-1, WaveList("*",";","DIMS:3")) Prompt source, "3D wave to show", popup, WaveList("*",";","DIMS:3") DoPrompt "pick wave", source if (V_flag) return 1 endif endif if (exists(source)!=1 || DimSize($source,2)<1) if (strlen(source)<1) print "no 3-d waves here" else printf "wave '%s' is un-usable\r",source endif return 1 endif Display /W=(86,104,456,430)/K=1 TextBox/N=text1/F=0/S=3/A=LB/X=-31/Y=-12 source String slice = UniqueName("slice_temp",1,0) String binName = GetDataFolder(3)+UniqueName("slice_binning",3,0) // variable to hold binning Variable nx=DimSize($source,0) Variable ny=DimSize($source,1) Variable nz=DimSize($source,2) Variable z0=DimOffset($source,2) Variable dz=DimDelta($source,2) Make/N=(nx,ny)/O/U/I $slice Note $slice note($source) if (SetXYtoANGER($slice)) // not an Anger detector, so copy the scales. SetScale/P x DimOffset($source,0),DimDelta($source,0),WaveUnits($source,0), $slice SetScale/P y DimOffset($source,1),DimDelta($source,1),WaveUnits($source,1), $slice endif Variable/G $binName=round(nz/12) NVAR binning=$binName binning = (nz<=25) ? 1 : binning AppendImage $slice ModifyImage $slice ctab= {*,*,Terrain,1} ModifyGraph margin(right)=7,gfSize=14,width={Aspect,1} ModifyGraph tick=2,mirror=1,minor=1,standoff=0 ModifyGraph lblMargin(left)=63, lblLatPos(left)=0 ModifyGraph axOffset(left)=4.86,axOffset(bottom)=-0.75 TextBox/N=text0/F=0/S=3/A=LB/X=-32/Y=-1.5 "\\Z160" Slider sliceSlider,pos={9,5},size={46,199},proc=upDateSliceProc Slider sliceSlider,limits={z0,(nz-1)*dz+z0,dz},value= z0,live= 0,tkLblRot= 90 Button down,pos={12,207},size={17,20},proc=PushButtonProc,title="<" Button up,pos={36,207},size={17,20},proc=PushButtonProc,title=">" SetVariable setBin,pos={3,232},size={61,15},title="bin",format="%d" SetVariable setBin,limits={1,nz,1},value= $binName SetWindow kwTopWin, hook=onWindowClose SetWindow kwTopWin, note="sliceWindow=1;" SetWindow kwTopWin, note+="source="+source+";" SetWindow kwTopWin, note+="sourceFull="+GetDataFolder(3)+source+";" SetWindow kwTopWin, note+="binName="+binName+";" SetWindow kwTopWin, note+="value="+num2str(z0)+";" SetWindow kwTopWin, note+="killVariable="+binName+";" SetWindow kwTopWin, note+="killWave="+GetDataFolder(3)+slice+";" DoWindow /T $WinName(0,1), WinName(0,1)+":"+source upDateSliceProc("",z0,1) End //Function ShowSliceGraph(source) // String source // if (exists(source)!=1) // source = "" // elseif (DimSize($source,2)<1) // source = "" // endif // if (strlen(source)<1) // source = StringFromList(ItemsInList(WaveList("*",";","DIMS:3"))-1, WaveList("*",";","DIMS:3")) // Prompt source, "3D wave to show", popup, WaveList("*",";","DIMS:3") // DoPrompt "pick wave", source // if (V_flag) // return 1 // endif // endif // if (exists(source)!=1 || DimSize($source,2)<1) // if (strlen(source)<1) // print "no 3-d waves here" // else // printf "wave '%s' is un-usable\r",source // endif // return 1 // endif // Display /W=(86,104,456,430)/K=1 // TextBox/N=text1/F=0/S=3/A=RB/X=-12.00/Y=-13.00 source // String slice = UniqueName("slice_temp",1,0) // String binName = GetDataFolder(3)+UniqueName("slice_binning",3,0) // variable to hold binning // Variable nx=DimSize($source,0) // Variable ny=DimSize($source,1) // Variable nz=DimSize($source,2) // Variable z0=DimOffset($source,2) // Variable dz=DimDelta($source,2) // Make/N=(nx,ny)/O/U/I $slice // Note $slice note($source) // SetXYtoANGER($slice) // Variable/G $binName=round(nz/12) // NVAR binning=$binName // binning = (nz<=25) ? 1 : binning // AppendImage $slice //// ModifyImage $slice ctab= {*,*,Rainbow,0} // ModifyImage $slice ctab= {*,*,Terrain,1} // ModifyGraph margin(right)=60,gfSize=14,width={Aspect,1} // ModifyGraph tick=2,mirror=1,minor=1,standoff=0 // ModifyGraph axOffset(left)=-4.1,axOffset(bottom)=-0.75 // TextBox/N=text0/F=0/S=3/A=RB/X=2.55/Y=3.64 "\\Z160" // Slider sliceSlider,pos={322,10},size={46,199},proc=upDateSliceProc // Slider sliceSlider,limits={z0,(nz-1)*dz+z0,dz},value= z0,live= 0,tkLblRot= 90 // Button down,pos={325,212},size={17,20},proc=PushButtonProc,title="<" // Button up,pos={349,212},size={17,20},proc=PushButtonProc,title=">" // SetVariable setBin,pos={316,237},size={61,15},title="bin",format="%d" // SetVariable setBin,limits={1,nz,1},value= $binName // SetWindow kwTopWin, hook=onWindowClose // SetWindow kwTopWin, note="sliceWindow=1;" // SetWindow kwTopWin, note+="source="+source+";" // SetWindow kwTopWin, note+="sourceFull="+GetDataFolder(3)+source+";" // SetWindow kwTopWin, note+="binName="+binName+";" // SetWindow kwTopWin, note+="value="+num2str(z0)+";" // SetWindow kwTopWin, note+="killVariable="+binName+";" // SetWindow kwTopWin, note+="killWave="+GetDataFolder(3)+slice+";" // upDateSliceProc("",z0,1) //End Function onWindowClose(infoStr) String infoStr if (stringmatch(StringByKey("EVENT",infoStr),"kill")) String graphName = StringByKey("WINDOW",infoStr) GetWindow $graphName, note String noteStr = S_value Wave image = ImageNameToWaveRef(graphName,StringFromList(0,ImageNameList(graphName,";"))) String wName = StringFromList(0,ImageNameList(graphName,";")) String binName = StringByKey("binName",noteStr,"=") RemoveImage $wName wName=StringByKey("killWave",noteStr,"=") KillWaves $wName String varName=StringByKey("killVariable",noteStr,"=") KillVariables $varName endif return 0 End Function upDateSliceProc(name,value,event) String name // name of this slider control Variable value // value of the slider (if value==-1 increase, value=-2 decrease) Variable event // bit field: 1=valueset, 2=mousedown, 4=mouse up,8=mouse moved if (!(event & 1)) return 0 endif String graphName = StringFromList(0,WinList("*",";","WIN:")) if (strlen(graphName)<1) DoAlert 0, "No Window for doing a slice update" return 1 endif GetWindow $graphName, note String noteStr = S_value String str, source Wave image = ImageNameToWaveRef(graphName,StringFromList(0,ImageNameList(graphName,";"))) source = StringByKey("sourceFull",noteStr,"=") Wave ws=$source if (!WaveExists(ws)) print "wave '"+source+"' not found" return 1 endif Variable nz=DimSize($source,2) Variable z0=DimOffset($source,2) Variable dz=DimDelta($source,2) Variable vmax = (nz-1)*dz+z0 NVAR binning = $StringByKey("binName",noteStr,"=") binning = max(1,min(binning,nz)) Variable upDownButton=0 // 0=normal, -1=down, +1=up if (stringmatch(name,"up")) upDownButton = 1 elseif (stringmatch(name,"down")) upDownButton = -1 endif if (upDownButton) Variable oldVal = NumberByKey("value",noteStr,"=") if (numtype(oldVal)) print "old value not found in note" return 1 endif value = oldVal + sign(upDownButton)*binning*dz if (valuevmax) beep endif endif value = max(z0,min(vmax,value)) Slider sliceSlider, value=value Variable ivalue = round((value-z0)/dz) Variable i,im=min(nz,ivalue+binning) image = 0 for (i=ivalue;i=vmax || value<=z0) return -1 endif return 0 End Function PushButtonProc(ctrlName) : ButtonControl String ctrlName if (!cmpstr(ctrlName,"up")) upDateSliceProc(ctrlName,1,1) endif if (!cmpstr(ctrlName,"down")) upDateSliceProc(ctrlName,-1,1) endif End Function/T sumUpWavelength(wName,range,thresh) String wName // name of wave to sum up String range // specifies channel ranges to sum up Variable thresh // only include pixels that are >= thresh, 0 includes all // if (strlen(range)<1) // range = "0-225,227-271" // endif Variable printIt=0 if (strlen(wName)<1 || strlen(range)<1 || numtype(thresh) || thresh<0) String wList = WaveList("*",";","DIMS:3") String fldrSav= GetDataFolder(1) SetDataFolder root:raw: wList += WaveList("*",";","DIMS:3") SetDataFolder fldrSav if (strlen(range)<1) sprintf range, "0-%d",DimSize($wName,2)-1 endif thresh = numtype(thresh) ? 0 : abs(thresh) Prompt wName, "name of wave to sum", popup, wList Prompt range, "range of channels (or wavelengths) to sum" Prompt thresh, "threshold for including a pixel, 0 includes all" DoPrompt "pick a wave", wName, range, thresh if (V_flag) return "" endif printIt = 1 printf " sumUpWavelength(\"%s\",\"%s\",%d)\r",wName,range,round(thresh) endif thresh = round(thresh) if (strsearch(range,".",0)>=0) // if a decimal point found, assume input is wavelength, convert printf "converted lambda range = '%s'", range range = LambdaRange2Channels(range) printf " to channel range '%s'\r", range endif String outName = wName+"sum" if (thresh>0) outName += "_"+num2istr(thresh) endif if (exists(wName)!=1) wName = "root:raw:"+wName endif Variable nx=DimSize($wName,0) Variable ny=DimSize($wName,1) Variable nz=DimSize($wName,2) String command, switches=makeWaveSwitches (NumberByKey("NUMTYPE", WaveInfo($wName,0))) sprintf command, "Make/O/N=(%d,%d)%s %s",nx,ny,switches,outName Execute command Wave out=$outName Wave in=$wName SetXYtoANGER(out) String noteStr = note(in) noteStr = ReplaceStringByKey("sumSource",noteStr,GetWavesDataFolder(in,2),"=") noteStr = ReplaceNumberByKey("threshSum",noteStr,thresh,"=") Note /K out Note out, noteStr range = expandRange(range,";") Duplicate /O in in_temp_ in_temp_ = (in_temp_[p][q][r]>=thresh) ? in_temp_[p][q][r] : 0 // zero all pixels below thresh out = 0 Variable k,m for (k=0;k64) // skip the first few lines starting with '#' skip -= 1 Close fid // the lines in the file contain the following: 'time delay (祍)' followed by (xSize*ySize) pixel values // there is one line for each time slice (zSize of them) // and the final line is end time of the last slice Variable start = ticks LoadWave/J/L={0,skip,0,0,0 }/M/V={"\t "," ",0, 0 } fileName Variable readTimeSec = (ticks-start)/60.15 Variable i = ItemsInList(S_waveNames) if (i!=1) Abort "the LoadWave command loaded "+num2str(i)+" waves" endif String wave0Name = StringFromList(0,S_waveNames) Wave wav = $wave0Name Variable N0,N1 N0 = DimSize(wav,0) N1 = DimSize(wav,1) Make/N=(N0)/O delays_temp_ delays_temp_ = wav[p][0] Redimension/n=(N0*N1) wav DeletePoints 0,N0, wav N1 -= 1 Redimension/n=(N0,N1) wav Variable Nx,Ny,Nz, wrongSize=0 // make sure image dimensions match globals if (xSize*ySize == N1) Nx = xSize Ny = ySize else Nx = round(sqrt(N1)) Ny = round(sqrt(N1)) wrongSize = 1 endif Nz = N0-1 wrongSize = wrongSize || (Nz !=zSize) if (wrongSize) Prompt Nx, "number of pixels along x direction, global="+num2str(xSize) Prompt Ny, "number of pixels along y direction, global="+num2str(ySize) Prompt Nz, "number of time slices in detector, global="+num2str(zSize) DoPrompt "Size of image in file does not match globals",Nx,Ny,Nz if (V_flag) DoAlert 0, "Leaving original values as they are" else xSize = Nx ySize = Ny zSize = Nz printf "updating global images sizes to xSize=%d, ySize=%d, zSize=%d\r",xSize,ySize,zSize endif endif Make/N=(Nx,Ny,Nz)/U/W/O $wName // unsigned 2-byte int Wave scd = $wName SetXYtoANGER(scd) scd = wav[r][Nx*q+p] // this does not copy the last line that is NaN's Note /K scd Note scd wNote printf "reading took %g sec, whole load of this file took %g sec\r",readTimeSec,(ticks-start)/60.15 KillWaves/Z $wave0Name UpdateGlobalsFromList(wNote) // update globals to reflect wave note values // check that delays matches the delays_temp_ Variable updateDelay = (exists("delays")!=1) // need to update (or make) delays if (!updateDelay) if (!compareWaves(delays,delays_temp_)) DoAlert 1, "delays do not match, update?" updateDelay = (V_flag==1) endif endif if (updateDelay) print "updating the 'delays' wave" Duplicate/O delays_temp_ delays SetScale d 0,0,"祍", delays endif KillWaves/Z delays_temp_ End Static Function/T ParametersInNote(x2mm,y2mm,xLeftmm,yLowermm,deta,detd,L1mm,Tzero,slitx,slitz,bx,slitWidth,surface) Variable x2mm // distance between pixels in x direction (mm) Variable y2mm // distance between pixels in y direction (mm) Variable xLeftmm // distance of first (left most) from center (mm) Variable yLowermm // distance of first (bottom most) from center (mm) Variable deta // angle from incident beam ot center of detector (clockwise deg.) Variable detd // distance from sample to detector center (mm) Variable L1mm // distance from sample to source (mm) Variable Tzero // time offset (祍ec), actual time is (time in file) -Tzero Variable slitx // x position of slit (mm) Variable slitz // z position of slit (mm) Variable bx // x-offset of neutron beam horizontally from correct axis (mm) Variable slitWidth // width of slit (mm) Variable surface // z-position of surface, when translation=0, surface is z-coordinate of the surface (mm) if (numtype(x2mm+y2mm+xLeftmm+yLowermm+deta+detd+L1mm+Tzero+slitx+slitz+slitWidth+surface)) // need to ask for values x2mm = numtype(x2mm) ? NumVarOrDefault("root:Packages:SCD:x2mm",1.52215) : x2mm y2mm = numtype(y2mm) ? NumVarOrDefault("root:Packages:SCD:y2mm",1.52535) : y2mm xLeftmm = numtype(xLeftmm) ? NumVarOrDefault("root:Packages:SCD:xLeftmm",-66.58399) : xLeftmm yLowermm = numtype(yLowermm) ? NumVarOrDefault("root:Packages:SCD:yLowermm",-85.97354) : yLowermm deta = numtype(deta) ? NumVarOrDefault("root:Packages:SCD:deta",-75) : deta detd = numtype(detd) ? NumVarOrDefault("root:Packages:SCD:detd",230) : detd L1mm = numtype(L1mm) ? NumVarOrDefault("root:Packages:SCD:L1mm",9371.1) : L1mm // 9378 Tzero = numtype(Tzero) ? NumVarOrDefault("root:Packages:SCD:Tzero",21.4) : Tzero slitx = numtype(slitx) ? NumVarOrDefault("root:Packages:SCD:slitx",-6.9) : slitx slitz = numtype(slitz) ? NumVarOrDefault("root:Packages:SCD:slitz",0) : slitz bx = numtype(bx) ? NumVarOrDefault("root:Packages:SCD:bx",0) : bx slitWidth = numtype(slitWidth) ? NumVarOrDefault("root:Packages:SCD:slitWidth",1) : slitWidth surface = numtype(surface) ? NumVarOrDefault("root:Packages:SCD:surface",0) : surface Prompt x2mm, "x distance between pixels (mm)" Prompt y2mm, "y distance between pixels (mm)" Prompt xLeftmm, "dist: left most pixel to center (mm)" Prompt yLowermm, " dist: lowest pixel to center (mm)" Prompt deta, "angle from incident beam to center of detector (clockwise degrees)" Prompt detd "sample to detector center distance (mm)" Prompt L1mm, "sample to source distance (mm)" Prompt Tzero, "time offset (祍ec), actual time is (time in file) -This" Prompt slitx, "x position of slit (mm)" Prompt slitz, "z position of slit (mm)" Prompt bx, "x-offset of neutron beam horizontally from correct axis (mm)" Prompt slitWidth, "width of slit (mm)" Prompt surface, "z-coord of surface (@ trans=0) (mm)" DoPrompt "Detector Size (mm not cm)",x2mm,x2mm,xLeftmm,yLowermm,deta,detd,L1mm,Tzero if (V_flag) DoAlert 0, "Default values have been used, if you don't like thisyou shouldn't have pushed cancel" endif DoPrompt "slit (mm not cm)",slitx,slitz,bx,slitWidth,surface if (V_flag) DoAlert 0, "Default values have been used, if you don't like thisyou shouldn't have pushed cancel" endif endif String list="" list = ReplaceNumberByKey("x2mm",list,x2mm,"=") list = ReplaceNumberByKey("y2mm",list,y2mm,"=") list = ReplaceNumberByKey("xLeftmm",list,xLeftmm,"=") list = ReplaceNumberByKey("yLowermm",list,yLowermm,"=") list = ReplaceNumberByKey("deta",list,deta,"=") list = ReplaceNumberByKey("detd",list,detd,"=") list = ReplaceNumberByKey("L1mm",list,L1mm,"=") list = ReplaceNumberByKey("Tzero",list,Tzero,"=") list = ReplaceNumberByKey("slitx",list,slitx,"=") list = ReplaceNumberByKey("slitz",list,slitz,"=") list = ReplaceNumberByKey("bx",list,bx,"=") list = ReplaceNumberByKey("slitWidth",list,slitWidth,"=") list = ReplaceNumberByKey("surface",list,surface,"=") return list End Static Function UpdateGlobalsFromList(values) String values // keyed list of values for resetting global values String var,varList = "x2mm;y2mm;xLeftmm;yLowermm;deta;detd;L1mm;Tzero;slitx;slitz;bx;slitWidth;surface" String wholeVar Variable val,i=-1 do i += 1 var = StringFromList(i,varList) if (strlen(var)<1) break endif val = NumberByKey(var,values,"=") if (numtype(val)) continue endif var = "root:Packages:SCD:"+var if (exists(var)!=2) continue endif NVAR global=$var global = val while(1) End //Static Function UpdateGlobalsFromList(values) // String values // keyed list of values for resetting global values // String var,varList = "x2mm;y2mm;xLeftmm;yLowermm;deta;detd;L1mm;Tzero;slitx;slitz;slitWidth;surface" // Variable val,i=0 // do // var = StringFromList(i,varList) // if (strlen(var)<1) // break // endif // var = "root:Packages:SCD:"+var // NVAR global=$(var+"root:Packages:SCD:") // if (!NVAR_Exists(global)) // continue // endif // val = NumberByKey(var,values,"=") // if (numtype(val)) // continue // endif // global = val // i += 1 // while(1) //End // ==================================================================== // ==================================================================== // =========================== End Reading IN ============================= // ==================================================================== // ==================================================================== // ==================================================================== // =========================== Start Panel Parameters ======================= // ==================================================================== Function DetectorGeometryWindow() if (strlen(WinList("DetectorGeometry", ";","WIN:64"))>1) DoWindow /F DetectorGeometry return 1 endif NewPanel/K=1/W=(23,49,459,367) DoWindow /C DetectorGeometry Variable deta_deg = NumVarOrDefault("root:Packages:SCD:deta",-90) Variable detd_mm = NumVarOrDefault("root:Packages:SCD:detd",250) Variable detWidth = NumVarOrDefault("root:Packages:SCD:xSize",100)*NumVarOrDefault("root:Packages:SCD:x2mm",1.5) // find current detector description Variable current = 1e20 String icurrent="set me" SVAR DetectorList = root:Packages:SCD:DetectorList String dectectorNames="",desc,item Variable i do item = StringFromList(i,DetectorList) desc = StringByKey("desc", item,"=",",") if (strlen(desc)>0) dectectorNames += desc+";" if (abs(NumberByKey("deta",item,"=",",")-deta_deg)+abs(NumberByKey("detd",item,"=",",")-detd_mm) < current) current = abs(NumberByKey("deta",item,"=",",")-deta_deg)+abs(NumberByKey("detd",item,"=",",")-detd_mm) if (current <.01) icurrent = desc endif endif endif i += 1 while(strlen(desc)>0) icurrent="set me" PopupMenu DetectorListPopup,pos={220,2},size={180,20},proc=PopMenuSetFromDetectorList,title="Presets" String cmd sprintf cmd, "PopupMenu DetectorListPopup,mode=1,popvalue=\"%s\",value=\"%s\"",icurrent,dectectorNames Execute cmd SetDrawLayer/K UserBack SetDrawLayer/K ProgBack RedrawDetectorGeometryInPanel(deta_deg,detd_mm,detWidth,1) SetDrawLayer ProgBack SetDrawEnv fsize= 16 DrawText 10,18,"Set Detector Geometry" SetDrawLayer UserBack SetVariable detaDisp,pos={22,28},size={86,17},title="deta",fSize=12,limits={-180,0,0},value=root:Packages:SCD:deta SetVariable detdDisp,pos={22,52},size={86,17},title="detd",fSize=12,format="%g",limits={0,500,0},value=root:Packages:SCD:detd SetVariable xSizeDisp,pos={149,28},size={86,17},title="xSize",fSize=12,format="%d",limits={0,4096,0},value= root:Packages:SCD:xSize SetVariable ySizeDisp,pos={147,52},size={86,17},title="ySize",fSize=12,format="%d",limits={0,4096,0},value=root:Packages:SCD:ySize SetVariable zSizeDisp,pos={147,76},size={86,17},title="zSize",fSize=12,format="%d",limits={0,4096,0},barmisc={0,1000},value=root:Packages:SCD:zSize SetVariable x2mmDisp,pos={259,28},size={118,17},title="x2mm",fSize=12,format="%g",limits={-5,5,0},value=root:Packages:SCD:x2mm SetVariable y2mmDisp,pos={261,52},size={118,17},title="y2mm",fSize=12,format="%g",limits={-5,5,0},value=root:Packages:SCD:y2mm SetVariable xLeftmmDisp,pos={249,75},size={126,17},title="xLeftmm",fSize=12,format="%g",limits={-Inf,0,0},value=root:Packages:SCD:xLeftmm SetVariable yLowermmDisp,pos={239,99},size={136,17},title="yLowermm",fSize=12,format="%g",limits={-Inf,0,0},value=root:Packages:SCD:yLowermm SetVariable L1mmDisp,pos={298,217},size={118,17},title="L1mm",fSize=12,format="%g",limits={10,Inf,0},value= root:Packages:SCD:L1mm SetVariable TzeroDisp,pos={300,241},size={118,17},title="Tzero",fSize=12,format="%g",limits={-Inf,Inf,0},value=root:Packages:SCD:Tzero SetVariable slitxDisp,pos={22,262},size={118,17},title="slitx",fSize=12,format="%g",limits={-Inf,0,0},value=root:Packages:SCD:slitx SetVariable slitzDisp,pos={22,286},size={118,17},title="slitz",fSize=12,format="%g",limits={-Inf,Inf,0},value=root:Packages:SCD:slitz SetVariable bxDisp,pos={170,262},size={78,17},title="bx",fSize=12,format="%g",limits={-Inf,Inf,0},value=root:Packages:SCD:bx SetVariable slitWidthDisp,pos={171,286},size={103,17},title="slitWidth",fSize=12,format="%g",limits={0,Inf,0},value=root:Packages:SCD:slitWidth SetVariable surfaceDisp,pos={294,286},size={103,17},title="surface",fSize=12,format="%g",limits={-Inf,Inf,0},value=root:Packages:SCD:surface return 0 End Function RedrawDetectorGeometryInPanel(deta,detd,detWidth,coords) Variable deta // detector angle (degrees) (=-75) Variable detd // detector distance (mm) ( =230) Variable detWidth // full width of detector (mm) ( =150) Variable coords // if true, draw coordinates if (strlen(WinList("DetectorGeometry", ";","WIN:64"))<2) return 1 elseif (numtype(deta+detd+detWidth+coords)) DoAlert 0,"invalid values in 'RedrawDetectorGeometryInPanel'" return 1 endif DoWindow /F DetectorGeometry Variable xc=130,yc=250 // origin is at top left Variable dist = detd/2 // distance of detector in screen units Variable detHalf=detWidth/4 // half width of detector in screen units SetDrawLayer/K UserBack SetDrawLayer UserBack Variable cosdeta=cos(deta*PI/180), sindeta=sin(deta*PI/180) Variable detx = xc-dist*cosdeta Variable dety = yc+dist*sindeta SetDrawEnv arrow=1 ; DrawLine xc,yc,detx,dety // arrow to detector center SetDrawEnv linethick=5 // This is the detector DrawLine detx+detHalf*sindeta,dety+detHalf*cosdeta,detx-detHalf*sindeta,dety-detHalf*cosdeta SetDrawEnv textxjust= (deta>-90) ? 2 : 0 DrawText detx,dety,"detector" DrawText xc+(detx-xc)/2+4, yc+(dety-yc)/2+2,"detd="+num2str(detd)+"mm" // make a polygon of an arc from 0 to -deta, radius is 1/2 detector distance Variable angle,rad=dist/2,N=10,i String arc sprintf arc, "DrawPoly %g,%g,1,1,{",xc-rad,yc for (i=0;i<=N;i+=1) angle = i*deta/N*PI/180 arc += num2str(xc-rad*cos(angle))+","+num2str(yc+rad*sin(angle)) if (i0) if (strlen(item)<5) return 1 endif NVAR deta = root:Packages:SCD:deta NVAR detd = root:Packages:SCD:detd NVAR L1 = root:Packages:SCD:L1mm NVAR Tzero = root:Packages:SCD:Tzero NVAR xSize = root:Packages:SCD:xSize NVAR ySize = root:Packages:SCD:ySize NVAR zSize = root:Packages:SCD:zSize NVAR x2mm = root:Packages:SCD:x2mm NVAR y2mm = root:Packages:SCD:y2mm NVAR xLeft = root:Packages:SCD:xLeftmm NVAR yLower = root:Packages:SCD:yLowermm deta = NumberByKey("deta",item,"=",",") detd = NumberByKey("detd",item,"=",",") detWidth = NumberByKey("x2mm",item,"=",",")*NumberByKey("xSize",item,"=",",") RedrawDetectorGeometryInPanel(deta,detd,detWidth,0) L1 = NumberByKey("L1",item,"=",",") Tzero = NumberByKey("Tzero",item,"=",",") xSize = NumberByKey("xSize",item,"=",",") ySize = NumberByKey("ySize",item,"=",",") zSize = NumberByKey("zSize",item,"=",",") x2mm = NumberByKey("x2mm",item,"=",",") y2mm = NumberByKey("y2mm",item,"=",",") xLeft = NumberByKey("xLeft",item,"=",",") yLower = NumberByKey("yLower",item,"=",",") End Function/T FindMatchToDetectorList(wav,DetectorList) Wave wav String DetectorList // SVAR DetectorList = root:Packages:SCD:DetectorList String noteStr = note(wav) Variable deta = NumberByKey("deta",noteStr,"=") Variable detd = NumberByKey("detd",noteStr,"=") Variable L1 = NumberByKey("L1mm",noteStr,"=") Variable Tzero = NumberByKey("Tzero",noteStr,"=") // Variable xSize = NumberByKey("xSize",noteStr,"=") // Variable ySize = NumberByKey("ySize",noteStr,"=") // Variable zSize = NumberByKey("zSize",noteStr,"=") Variable x2mm = NumberByKey("x2mm",noteStr,"=") Variable y2mm = NumberByKey("y2mm",noteStr,"=") Variable xLeft = NumberByKey("xLeftmm",noteStr,"=") Variable yLower = NumberByKey("yLowermm",noteStr,"=") // Variable slitx = NumberByKey("slitx",noteStr,"=") // Variable slitz = NumberByKey("slitz",noteStr,"=") // Variable bx = NumberByKey("bx",noteStr,"=") // Variable slitWidth = NumberByKey("slitWidth",noteStr,"=") Variable diff, i String item for (i=0;i=48 && nchar<=57 || nchar==46) // a digit, or '.' num += char // part of a number, accumulate else if (strlen(num)>0) // done accumulating number, so include it if (dashOn) Crange += num2str(ceil(lambda2chan(str2num(num)))) else Crange += num2str(floor(lambda2chan(str2num(num)))) endif dashOn = (nchar==45) num = "" endif Crange += char endif endfor if (strlen(num)>0) // include any trailing numbers Crange += num2istr(lambda2chan(str2num(num))) num = "" endif return Crange End Function chan2lambda(index) // convert time channel to wavelength () Variable index // index into the time part of the data Variable printIt=0 NVAR detd = root:Packages:SCD:detd // distance from sample to detector center (mm) NVAR L1mm = root:Packages:SCD:L1mm // distance from sample to source (mm) NVAR Tzero = root:Packages:SCD:Tzero // time offset (祍ec), actual time is (time in file) -Tzero NVAR zSize = root:Packages:SCD:zSize // number of time slices if (numtype(index)) index = 0 Prompt index, "channel number (range 0-"+num2istr(zSize)+")" DoPrompt "pick a channel", index if (V_flag) return NaN endif printIt = 1 endif Variable distance = (L1mm + detd)/1000 // total path length of neutrons (m) Variable t0=Tzero // delay in 祍ec Variable lam1, lam2, lam // wavelengths at either end of the channel Wave delays = delays // prefer local values if (!WaveExists(delays)) Wave delays = root:raw:delays endif if (!WaveExists(delays)) Abort "cannot find the 'delays' wave in chan2lambda("+num2str(index)+")" endif if (numtype(index)) return NaN endif index = max(index,0) index = min(index,DimSize(delays,0)-2) lam1 = velocity2lambda(distance / (delays[index]+t0) * 1e6) lam2 = velocity2lambda(distance / (delays[index+1] +t0)* 1e6) lam = (lam1+lam2)/2 if (printIt) printf "for channel %d, the wavelength is %g ()\r", index, lam endif return lam End Function lambda2chan(lambda) // convert wavelength () to time channel, returns NaN on error Variable lambda // wavelength () Variable printIt=0 if (numtype(lambda) || lambda<0) lambda = 1.8 Prompt lambda, "wavelength ()" DoPrompt "pick a wavelength", lambda if (V_flag) return NaN endif printIt = 1 endif NVAR detd = root:Packages:SCD:detd // distance from sample to detector center (mm) NVAR L1mm = root:Packages:SCD:L1mm // distance from sample to source (mm) NVAR Tzero = root:Packages:SCD:Tzero // time offset (祍ec), actual time is (time in file) -Tzero // Variable delay = distance / lambda2velocity(lambda) // neutron velocity from the wavelength Variable distance = (L1mm + detd)/1000 // path length of neutrons (m) Variable delay = distance/lambda2velocity(lambda) // neutron velocity from the wavelength delay = delay*1e6 - Tzero // true delay (祍ec) Wave delays = delays // prefer local values if (!WaveExists(delays)) Wave delays = root:raw:delays endif Variable index = BinarySearchInterp(delays,delay)// index into the time part of the data if (!printIt) return index // just return value, no matter what elseif (delaydelays[inf]) print "wavelength is too long" DoAlert 0, "wavelength is too long" elseif (numtype(index)) print "something is very wrong computing wavelength" Abort "something is very wrong computing wavelength" elseif (printIt) printf "for a wavelength of %g (), the nearest channel is %.1f\r", lambda,index endif return index End Static Function/T WaveList_NxN(n0,n1) Variable n0,n1 String wList=WaveList("*",";","DIMS:2"), NxNList="", item Variable i for (i=0;i2.1) // common factors of 3 or greater, definitely not lowest order KillWaves/Z h_temp_ return 0 elseif (factor <1.1) // no common factors KillWaves/Z h_temp_ return 1 endif // at this point, it is allowed, and factor is 2, so check if factor of 2 is needed h_temp_ /= 2 Variable i = !allowBCC(h_temp_) KillWaves/Z h_temp_ return i End // Determine the angle between two sets of G vectors whose diffracted spots are at cursors A and B // Note, that this is not the angle between the cursors themselves, but between the G vectors // Note that +z is along beam, +y is up, and +x is horizontally perp to beam (points to where you stand). Function angleBetweenGofCursors() // determine angle of G vector between cursor A and B Variable x0 = hcsr(A) // pixel poxitions on detector (mm from center) Variable y0 = vcsr(A) // not beamline coords Variable x1 = hcsr(B) Variable y1 = vcsr(B) String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";")))) Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta") Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd") Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx") Variable xp, zp // pixel position (in horizontal plane) Make/O/D/N=3 kihat={0,0,1}, kfhat, G0,G1 // kihat points downstream, from the sampe away from the source pixelXZfromDetector(deta,detd,x0,xp,zp) // calculate pixel coordinates xp,zp kfhat = {xp-bx,y0,zp} // vertical height (y-direction) is y0 normalize(kfhat) G0 = (kfhat-kihat) // points half way between kf and -ki pixelXZfromDetector(deta,detd,x1,xp,zp) // calculate pixel coordinates xp,zp kfhat = {xp-bx,y1,zp} // vertical height (z-direction) is y1 normalize(kfhat) G1 = (kfhat-kihat) Variable angle = angleBetweenVectors(G0,G1)*180/PI KillWaves G0,G1, kihat,kfhat return angle End //Function angleBetweenGofCursors() // determine angle of G vector between cursor A and B // Variable x0 = hcsr(A) // pixel poxitions on detector (mm from center) // Variable y0 = vcsr(A) // not beamline coords // Variable x1 = hcsr(B) // Variable y1 = vcsr(B) // // String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";")))) // Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta") // Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd") // Variable xp, zp // pixel position (in horizontal plane) // Make/O/D/N=3 kihat={0,0,1}, kfhat, G0,G1 // // kihat points downstream, from the sampe away from the source // pixelXZfromDetector(deta,detd,x0,xp,zp) // calculate pixel coordinates xp,zp // kfhat = {xp,y0,zp} // vertical height (y-direction) is y0 // normalize(kfhat) // G0 = (kfhat-kihat) // points half way between kf and -ki // // pixelXZfromDetector(deta,detd,x1,xp,zp) // calculate pixel coordinates xp,zp // kfhat = {xp,y1,zp} // vertical height (z-direction) is y1 // normalize(kfhat) // G1 = (kfhat-kihat) // // Variable angle = angleBetweenVectors(G0,G1)*180/PI // KillWaves G0,G1, kihat,kfhat // return angle //End Function pixelXZfromDetector(deta,detd,xDet,xp,zp) Variable deta // angle from exit beam to detector center (measured counter clock wise) (deg) Variable detd // distance from sample to detector center (mm) Variable xDet // horizontal distance along the detector from detector center to pixel (mm) Variable &xp, &zp // pixel position in beamline coordinates (in horizontal plane) (mm) Variable cosDeta = cos(deta*PI/180) Variable sinDeta = sin(deta*PI/180) xp = detd*sinDeta - xDet*cosDeta // pixel position (in horizontal plane) zp = detd*cosDeta + xDet*sinDeta return NaN End Function RemoveCommonFactors(vec) // if elements of vec are integres, remove all common factors Wave vec Variable N=numpnts(vec) WaveStats/Q vec V_max = abs(V_max) if (V_max>2000) Abort "elements of vector way too big for factoring" endif Variable i for (i=V_max; i>1; i-=1) if (mod(vec[0],i)==0 && mod(vec[1],i)==0 && mod(vec[2],i)==0) vec /= i endif endfor Variable vm=V_max WaveStats/Q vec return vm/V_max End //Function MatrixRotate(index,angle,mat) // rotates matrix about x,y,z for index = 0,1,2 // Variable index // Variable angle // (radians) // Wave mat // // Variable N = DimSize(mat,0) // if (DimSize(mat,1)!=N && DimSize(mat,2)!=0) // Abort "can only rotate square 2-d matricies" // endif // // Variable cosa = cos(angle) // Variable sina = sin(angle) // // Duplicate/O mat rot_temp_ // rot_temp_ = 0 // rot_temp_[index][index] = 1 // // Variable i0 = mod(index+1,N) // Variable i1 = mod(index+2,N) // rot_temp_[i0][i0] = cosa // rot_temp_[i1][i1] = cosa // rot_temp_[i0][i1] = sina // rot_temp_[i1][i0] = -sina // // MatrixMultiply mat, rot_temp_ // Wave M_product=M_product // mat = M_product // KillWaves/Z rot_temp_ //End // ==================================================================== // ============================ End Utility ============================== // ==================================================================== // ==================================================================== // ==================================================================== // ==================================================================== // =========================== Start Initialize ============================= // ==================================================================== Window Panel_SCD_geometry() : Panel PauseUpdate; Silent 1 // building window... if (strlen(WinList("SCDgeometry", ";","WIN:64"))>1) DoWindow /F SCDgeometry return endif NewPanel /K=1 /W=(451,52,760,198) DoWindow /C SCDgeometry SetDrawLayer UserBack SetDrawEnv fillfgc= (1,52428,26586) DrawOval 112,121,123,132 SetDrawEnv linethick= 2,arrow= 2 DrawLine 123,126,275,126 DrawText 154,123,"source to sample is L1" SetDrawEnv linethick= 5 DrawLine 121,7,19,62 SetDrawEnv arrow= 1 DrawLine 115,122,69,36 SetDrawEnv dash= 3 DrawLine 110,126,10,126 DrawText 90,73,"detd" DrawText 30,106,"裠eta" SetDrawEnv arrow= 1,fillpat= 0 DrawPoly 59,127,1,1,{85,148,89,130,97,116,111,104,120,100,120,100} DrawText 43,25,"detector" SetDrawEnv linethick= 2,arrow= 1 DrawLine 233,40,233,87 SetDrawEnv linethick= 2,arrow= 1 DrawLine 234,40,185,40 DrawText 213,87,"+x" DrawText 171,46,"+z" DrawText 246,44,"+y" SetDrawEnv linethick= 2,fillpat= 0,fillfgc= (1,52428,26586) DrawOval 232,31,243,42 SetDrawEnv fillfgc= (0,0,0) DrawOval 235,35,240,39 EndMacro Function SCDinitPackage(doDialog) Variable doDialog // if true, do the dialog even if all values look OK // this routine sets up global variables that describe the size and position of the detector if (exists("NeutronDataInitPackage")!=6 || exists("root:Packages:Neutron:Isotope")!=1 || exists("root:Packages:Neutron:Abs_xs")!=1) Execute/P "INSERTINCLUDE \"Neutron\", version>=1.0" Execute/P "COMPILEPROCEDURES" NeutronDataInitPackage() endif // fwhm is only used in the reolution stuff // Execute "FWHM",Execute/P "INSERTINCLUDE \"FWHM\"";Execute/P "COMPILEPROCEDURES" NewDataFolder /O root:Packages NewDataFolder /O root:Packages:SCD // Fri Sep 12, 2003 10:01:35 AM // Jon, // Here are the current calibration parameters. // // DN DETA DETD L1 Tzero x2cm y2cm xLeft yLower Description // 17 -75.00 23.00 937.11 21.4 0.152215 0.152535 -6.658399 -8.597354 5-SCDcal 2003-08-05 // 19 -120.00 18.00 937.11 21.4 0.145137 0.144757 -6.551019 -8.185886 5-SCDcal 2003-08-05 // // Where DETA the detector angle in degrees, DETD is the sample-to-detector // distance in centimeters, L1 is the sample-to-moderator distance in // centimeters, x2cm and y2cm are channel widths in cm, xLeft and yLower // are the distance of the left and lower edges from the "center" of the // detector, which has coordinates of 0.0,0.0 in centimeter units. // // the coordinate of the [0,0] pixel is (xLeft,yLower) // // Art // IPNS, Bldg. 360 Ph: 630-252-3465 // Much OLDER values // // from instprm.dat // phi=75,chi=90,omega=45 // // DN DETA DETD L1 Tzero x2cm y2cm xLeft yLower // 17 -90.00 25.00 933.68 0.0 0.1 0.1 -7.5 -7.5 if (exists("root:Packages:SCD:DetectorList")!=2) String/G root:Packages:SCD:DetectorList = "" // list of detector parameters SVAR DetectorList=root:Packages:SCD:DetectorList DetectorList = AddListItem("deta=-90,detd=320,L1=9437.1,Tzero=5.5,xSize=85,ySize=85,zSize=120,x2mm=2.9267,y2mm=3.0238,xLeft=-126.157,yLower=-124.4681,desc=first single 90 detector",DetectorList) DetectorList = AddListItem("deta=-120,detd=180,L1=9378,Tzero=0,xSize=100,ySize=100,zSize=272,x2mm=1.5,y2mm=1.5,xLeft=-75,yLower=-75,desc=120 nominal detector",DetectorList) DetectorList = AddListItem("deta=-75,detd=230,L1=9378,Tzero=0,xSize=100,ySize=100,zSize=272,x2mm=1.5,y2mm=1.5,xLeft=-75,yLower=-75,desc=75 nominal detector",DetectorList) DetectorList = AddListItem("deta=-120,detd=180,L1=9371.1,Tzero=21.4,xSize=100,ySize=100,zSize=272,x2mm=1.45137,y2mm=1.44757,xLeft=-65.51019,yLower=-81.85886,desc=120 2003-08-05",DetectorList) DetectorList = AddListItem("deta=-75,detd=230,L1=9371.1,Tzero=21.4,xSize=100,ySize=100,zSize=272,x2mm=1.52215,y2mm=1.52535,xLeft=-66.58399,yLower=-85.97354,desc=75 2003-08-05",DetectorList) endif if (exists("root:Packages:SCD:xSize")!=2) Variable/G root:Packages:SCD:xSize = NaN // x dimension of detector (pixels) endif if (exists("root:Packages:SCD:ySize")!=2) Variable/G root:Packages:SCD:ySize = NaN // y dimension of detector (pixels) endif if (exists("root:Packages:SCD:zSize")!=2) Variable/G root:Packages:SCD:zSize = NaN // number of time slices endif if (exists("root:Packages:SCD:x2mm")!=2) Variable/G root:Packages:SCD:x2mm = NaN // distance between pixels in x direction (mm) endif if (exists("root:Packages:SCD:y2mm")!=2) Variable/G root:Packages:SCD:y2mm = NaN // distance between pixels in y direction (mm) endif if (exists("root:Packages:SCD:xLeftmm")!=2) Variable/G root:Packages:SCD:xLeftmm = NaN // distance of first (left most) from center (mm) endif if (exists("root:Packages:SCD:yLowermm")!=2) Variable/G root:Packages:SCD:yLowermm = NaN // distance of first (bottom most) from center (mm) endif if (exists("root:Packages:SCD:phi")!=2) Variable/G root:Packages:SCD:deta = NaN // angle from incident beam to center of detector (clockwise deg.) endif if (exists("root:Packages:SCD:detd")!=2) Variable/G root:Packages:SCD:detd = NaN // distance from sample to detector center (mm) endif if (exists("root:Packages:SCD:L1mm")!=2) Variable/G root:Packages:SCD:L1mm = NaN // distance from sample to source (mm) endif if (exists("root:Packages:SCD:Tzero")!=2) Variable/G root:Packages:SCD:Tzero = NaN // time offset (祍ec), actual time is (time in file) -Tzero endif if (exists("root:Packages:SCD:slitx")!=2) Variable/G root:Packages:SCD:slitx = NaN // x position of slit (mm) endif if (exists("root:Packages:SCD:slitz")!=2) Variable/G root:Packages:SCD:slitz = NaN // z position of slit (mm) endif if (exists("root:Packages:SCD:bx")!=2) Variable/G root:Packages:SCD:bx = NaN // x-offset of neutron beam horizontally from correct axis (mm) endif if (exists("root:Packages:SCD:slitWidth")!=2) Variable/G root:Packages:SCD:slitWidth = NaN // width of slit (mm) endif if (exists("root:Packages:SCD:surface")!=2) Variable/G root:Packages:SCD:surface = NaN // z-coord of surface (when trans=0) (mm) endif if (exists("root:Packages:SCD:rangeList")!=2) String/G root:Packages:SCD:rangeList = "8461-8491;8492-8522;8523-8538;" endif if (exists("root:Packages:SCD:materialsList")!=2) String/G root:Packages:SCD:materialsList = "Cu;Nb;" endif if (exists("root:Packages:SCD:aoList")!=2) String/G root:Packages:SCD:aoList = "" endif if (exists("root:Packages:SCD:BravaisList")!=2) String/G root:Packages:SCD:BravaisList = "" endif NVAR xSize = root:Packages:SCD:xSize // x dimension of detector (pixels) NVAR ySize = root:Packages:SCD:ySize // y dimension of detector (pixels) NVAR zSize = root:Packages:SCD:zSize // number of time slices NVAR x2mm = root:Packages:SCD:x2mm // distance between pixels in x direction (mm) NVAR y2mm = root:Packages:SCD:y2mm // distance between pixels in y direction (mm) NVAR xLeftmm = root:Packages:SCD:xLeftmm // distance of first (left most) from center (mm) NVAR yLowermm = root:Packages:SCD:yLowermm// distance of first (bottom most) from center (mm) NVAR deta = root:Packages:SCD:deta // angle from incident beam ot center of detector (clockwise deg.) NVAR detd = root:Packages:SCD:detd // distance from sample to detector center (mm) NVAR L1mm = root:Packages:SCD:L1mm // distance from sample to source (mm) NVAR Tzero = root:Packages:SCD:Tzero // time offset (祍ec), actual time is (time in file) -Tzero NVAR slitx = root:Packages:SCD:slitx // x position of slit (mm) NVAR slitz = root:Packages:SCD:slitz // z position of slit (mm) NVAR bx = root:Packages:SCD:bx // x position of neutron beam (beam not on z-axis) (mm) NVAR slitWidth = root:Packages:SCD:slitWidth// width of slit (mm) NVAR surface = root:Packages:SCD:surface // z-coord of surface (when trans=0) (mm) SVAR materialsList = root:Packages:SCD:materialsList // list of all materials in the sample doDialog += numtype(xSize)|| numtype(ySize)|| numtype(zSize) doDialog += numtype(x2mm)|| numtype(y2mm)|| numtype(xLeftmm)|| numtype(yLowermm) doDialog += numtype(deta)|| numtype(detd)|| numtype(L1mm)|| numtype(Tzero) doDialog += numtype(slitx+slitz+slitWidth+bx) doDialog += (strlen(materialsList)<1) if (doDialog) xSize = numtype(xSize) ? 100 : xSize ySize = numtype(ySize) ? 100 : ySize zSize = numtype(zSize) ? 272 : zSize x2mm = numtype(x2mm) ? 1.52215 : x2mm y2mm = numtype(y2mm) ? 1.52535 : y2mm xLeftmm = numtype(xLeftmm) ? -66.58399 : xLeftmm yLowermm = numtype(yLowermm) ? -85.97354 : yLowermm deta = numtype(deta) ? -75 : deta detd = numtype(detd) ? 230 : detd L1mm = numtype(L1mm) ? 9371.1 : L1mm // 9336.8 Tzero = numtype(Tzero) ? 21.4 : Tzero slitx = numtype(slitx) ? -6.9 : slitx slitz = numtype(slitz) ? 0.00 : slitz bx = numtype(bx) ? 0.0 : bx slitWidth = numtype(slitWidth) ? 1.00 : slitWidth surface = numtype(surface) ? 0 : surface Variable x_Size=xSize, y_Size=ySize, z_Size=zSize Variable dxpixel_mm=x2mm, dypixel_mm=y2mm Variable xLeft_mm=xLeftmm, yLower_mm=yLowermm Variable deta_deg=deta, detd_mm=detd, L1_mm=L1mm, Tzero_usec=Tzero Variable slitx_mm=slitx, slitz_mm=slitz, bx_mm=bx, slitWidth_mm=slitWidth,surface_mm=surface String materials_List=materialsList Prompt x_Size, "number of pixels along x direction" Prompt y_Size, "number of pixels along y direction" Prompt z_Size, "number of time slices in detector" Prompt dxpixel_mm, "distance between pixels in x (mm)" Prompt dypixel_mm, "distance between pixels in y (mm)" Prompt xLeft_mm, "dist: left most pixel to center (mm)" Prompt yLower_mm, " dist: lowest pixel to center (mm)" Prompt deta_deg, "angle from incident beam to center of detector (clockwise degrees)" Prompt detd_mm, "distance from sample to detector center (mm)" Prompt L1_mm, "distance from sample to source (mm)" Prompt Tzero_usec, "time offset (祍ec), actual time is (time in file) -This" Prompt slitx_mm, "x position of slit (mm)" Prompt slitz_mm, "z position of slit (mm)" Prompt bx_mm, "offset of neutron beam in x-direcion from correct axis (mm)" Prompt slitWidth_mm, "width of slit (mm)" Prompt surface_mm, "z-coord of surface (@ trans=0) (mm)" Prompt materials_List, "list of all materials, semi-colon separated" DoPrompt "Detector Size (mm not cm)",x_Size,y_Size,z_Size,dxpixel_mm,dypixel_mm,xLeft_mm,yLower_mm if (V_flag) DoAlert 0, "Default values have been set, if you don't like this, rerun 'SCDinitPackage(1)'" endif xSize = x_Size ySize = y_Size zSize = z_Size x2mm = dxpixel_mm y2mm = dypixel_mm xLeftmm = xLeft_mm yLowermm = yLower_mm if (V_flag) return 1 endif DoPrompt "Detector Placement (mm not cm)",deta_deg,detd_mm,L1_mm,Tzero_usec if (V_flag) DoAlert 0, "Default values have been set, if you don't like this, rerun 'SCDinitPackage(1)'" endif deta = deta_deg detd = detd_mm L1mm = L1_mm Tzero = Tzero_usec DoPrompt "Slit Placement (mm not cm)",slitx_mm,slitz_mm,bx_mm,slitWidth_mm,surface_mm if (V_flag) DoAlert 0, "Default values have been set, if you don't like this, rerun 'SCDinitPackage(1)'" endif slitx = slitx_mm slitz = slitz_mm bx = bx_mm slitWidth = slitWidth_mm surface = surface_mm DoPrompt "materials",materials_List if (V_flag) DoAlert 0, "Default values have been set, if you don't like this, rerun 'SCDinitPackage(1)'" endif materialsList = materials_List endif // now make sure that lattice constants and types are set properly SVAR aoList = root:Packages:SCD:aoList // lattice constants list of all materials in materialsList SVAR BravaisList = root:Packages:SCD:BravaisList // lattice type list of all materials in materialsList Variable i, N=ItemsInList(materialsList) if (!doDialog && N==ItemsInList(aoList) && N==ItemsInList(BravaisList)) // should be OK return 0 endif aoList="" BravaisList="" for (i=0;i