Changeset 1212
 Timestamp:
 Sep 18, 2019 2:31:55 PM (3 years ago)
 Location:
 sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS
 Files:

 3 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/AvgGraphics.ipf
r1208 r1212 32 32 //ok, create the globals, fill the keyword string with all possible values (default) 33 33 String/G root:myGlobals:Drawing:gDrawInfoStr = "AVTYPE=Circular;PHI=0;DPHI=0;WIDTH=0;SIDE=both;" 34 root:myGlobals:Drawing:gDrawInfoStr += "QCENTER=0;QDELTA=0; "34 root:myGlobals:Drawing:gDrawInfoStr += "QCENTER=0;QDELTA=0;RATIOAXES=1;" 35 35 Variable/G root:myGlobals:Drawing:gDrawPhi =0 36 36 Variable/G root:myGlobals:Drawing:gDrawWidth = 1 … … 133 133 break 134 134 case "Elliptical": 135 EllipticalAverageTo1D(type)136 break137 135 case "Circular": 138 136 case "Sector": 139 //circular or sector137 //circular, elliptical, or sector 140 138 CircularAverageTo1D(type) //graph is drawn here 141 139 break … … 309 307 SetVariable width_p,disable=yes 310 308 SetVariable RAxes_p,disable=no 311 popupmenu sides,disable= no309 popupmenu sides,disable=yes 312 310 break 313 311 default: // optional default expression executed … … 350 348 SetVariable DPhi_p,help={"Enter the +/ range (0,45) of azimuthal angles to be included in the average. The bounding angles will be drawin in blue."} 351 349 SetVariable DPhi_p,limits={0,90,1},value= root:myGlobals:Drawing:gDrawDPhi 352 SetVariable RAxes_p,pos={166,176},size={110,15},proc= DeltaPhiSetVarProc,title="Axis Ratio"350 SetVariable RAxes_p,pos={166,176},size={110,15},proc=AxisRatioSetVarProc,title="Axis Ratio" 353 351 SetVariable RAxes_p,help={"Enter the ratio of the minor to major axes for an isointensity contour. By definition, this value should be less than 1."} 354 352 SetVariable RAxes_p,limits={0,1,0.001},value= root:myGlobals:Drawing:gDrawRAxes … … 491 489 variable RAxes = NumberByKey("RATIOAXES",drawStr,"=",";") 492 490 491 rr=0 492 gg=50000 493 bb=0 494 493 495 Make/O/N=(128,128) ellipsewave 494 496 variable b = 64 495 497 variable a = RAxes * b 496 DrawAnEllipsoid("ellipsewave",x0,y0,a,b,phi )498 DrawAnEllipsoid("ellipsewave",x0,y0,a,b,phi,rr,gg,bb) 497 499 498 500 Return 0 //exit the Draw routine … … 826 828 // OUTPUTS: 827 829 // make all changes directly in xywave 828 Function DrawAnEllipsoid(wname,xo,yo,a,b,phi )830 Function DrawAnEllipsoid(wname,xo,yo,a,b,phi,rr,gg,bb) 829 831 string wname 830 variable xo,yo,a,b,phi 832 variable xo,yo,a,b,phi,rr,gg,bb 831 833 int ii 832 834 833 835 wave xywave = $wname 834 836 835 make/ FREE/n=361 urx,ury836 make/ n=361 xw,yw = 0837 make/O/FREE/N=361 urx,ury 838 make/O/N=361 xw,yw = 0 837 839 838 840 variable cosalpha = cos(phi*Pi/180), sinalpha = sin(phi*Pi/180) … … 848 850 xywave[xw,yw] = 1 849 851 850 rr=0851 gg=50000852 bb=0853 854 852 //SetDrawEnv/W=SANS_Data xcoord= bottom,ycoord= left,linefgc= (rr,gg,bb),linethick= (thick),fillpat=0 855 853 //AppendImage/W=SANS_Data ellipsewave … … 944 942 ControlInfo/W=Average_Panel RAxes_p 945 943 Variable val = V_Value 946 SVAR tempStr = root:myGlobals:Drawing:gDraw RAxes944 SVAR tempStr = root:myGlobals:Drawing:gDrawInfoStr 947 945 String newStr = ReplaceNumberByKey("RATIOAXES", tempStr, val, "=", ";") 948 946 949 String/G root:myGlobals:Drawing:gDraw RAxes= newStr947 String/G root:myGlobals:Drawing:gDrawInfoStr = newStr 950 948 951 949 //redraw the angles 
sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/CircSectAve.ipf
r902 r1212 32 32 33 33 SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr //this is the list that has it all 34 Variable isCircular = 0 34 Variable isCircular = 0, isElliptical = 0 35 35 36 36 if( cmpstr("Circular",StringByKey("AVTYPE",keyListStr,"=",";")) ==0) 37 37 isCircular = 1 //set a switch for later 38 Endif 39 if( cmpstr("Elliptical",StringByKey("AVTYPE",keyListStr,"=",";")) ==0) 40 isElliptical = 1 //set a switch for later 38 41 Endif 39 42 … … 44 47 // 45 48 Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist,dr,ddr 46 Variable lambda,trans 49 Variable lambda,trans,raxes 47 50 WAVE reals = $(destPath + ":RealsRead") 48 51 WAVE/T textread = $(destPath + ":TextRead") … … 88 91 phi_y = sin(phi_rad) 89 92 Endif 93 if(isElliptical) 94 raxes = NumberByKey("RATIOAXES",keyListStr,"=",";") 95 EndIf 90 96 91 97 /// data wave is data in the current folder which was set at the top of the function … … 145 151 //BEGIN AVERAGE ********** 146 152 Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1 147 Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p 153 Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p,rho 148 154 149 155 // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too) … … 178 184 do 179 185 dyy = dy + (kk  fd)*sy/3 180 if(isCircular) 181 //circular average, use all pixels 182 //(increment) 183 nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2) 186 if(isCircular  isElliptical) 187 //use all pixels 188 if (isElliptical) 189 rho = atan(dyy/dxx)  phi_rad 190 nq = IncrementEllipticalPixel(data_pixel,ddr,dxx,dyy,rho,raxes, aveint,dsq,ncells,nq,nd2) 191 else 192 //circular average 193 nq = IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2) 194 EndIf 184 195 else 185 196 //a sector average  determine azimuthal angle … … 396 407 End 397 408 409 //returns nq, new number of qvalues 410 //arrays aveint,dsq,ncells are also changed by this function 411 // 412 Function IncrementEllipticalPixel(dataPixel,ddr,dxx,dyy,rho,raxes,aveint,dsq,ncells,nq,nd2) 413 Variable dataPixel,ddr,dxx,dyy,rho,raxes 414 Wave aveint,dsq,ncells 415 Variable nq,nd2 416 417 Variable irCircular,ir 418 419 irCircular = (sqrt(dxx*dxx+dyy*dyy)/ddr) 420 ir = irCircular*sqrt(cos(rho)*cos(rho) + raxes*raxes*sin(rho)*sin(rho))+1 421 if (ir>nq) 422 nq = ir //resets maximum number of qvalues 423 endif 424 aveint[ir1] += dataPixel/nd2 //ir1 must be used, since ir is physical 425 dsq[ir1] += dataPixel*dataPixel/nd2 426 ncells[ir1] += 1/nd2 427 428 Return nq 429 End 430 398 431 //function determines azimuthal angle dphi that a vector connecting 399 432 //center of detector to pixel makes with respect to vector 
sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/RectAnnulAvg.ipf
r1208 r1212 660 660 Return 0 661 661 End 662 663 //performs an elliptical average using isointensity contours, centered on a664 //specified qvalue (Intensity vs. angle)665 //the parameters in the global keywordstring must have already been set somewhere666 //either directly or from the protocol667 //668 //the input (data in the "type" folder) must be on linear scale  the calling routine is669 //responsible for this670 //averaged data is written to the data folder and plotted. data is not written671 //to disk from this routine.672 //673 Function EllipticalAverageTo1D(type)674 String type675 676 SVAR keyListStr = root:myGlobals:Protocols:gAvgInfoStr677 678 //type is the data type to do the averaging on, and will be set as the current folder679 //get the current displayed data (so the correct folder is used)680 String destPath = "root:Packages:NIST:"+type681 682 Variable xcenter,ycenter,x0,y0,sx,sx3,sy,sy3,dtsize,dtdist683 Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,xbm,ybm,ii684 Variable rc,delr,rlo,rhi,dphi,nphi,dr685 Variable lambda,trans686 Wave reals = $(destPath + ":RealsRead")687 688 // center of detector, for nonlinear corrections689 NVAR pixelsX = root:myGlobals:gNPixelsX690 NVAR pixelsY = root:myGlobals:gNPixelsY691 692 xcenter = pixelsX/2 + 0.5 // == 64.5 for 128x128 Ordela693 ycenter = pixelsY/2 + 0.5 // == 64.5 for 128x128 Ordela694 695 // beam center, in pixels696 x0 = reals[16]697 y0 = reals[17]698 //detector calibration constants699 sx = reals[10] //mm/pixel (x)700 sx3 = reals[11] //nonlinear coeff701 sy = reals[13] //mm/pixel (y)702 sy3 = reals[14] //nonlinear coeff703 704 dtsize = 10*reals[20] //det size in mm705 dtdist = 1000*reals[18] // det distance in mm706 lambda = reals[26]707 708 Variable phi = NumberByKey("PHI",keyListStr,"=",";")709 Variable ratio = NumberByKey("RATIOAXES",keyListStr,"=",";")710 711 dr = 1 //minimum annulus width, keep this fixed at one712 NVAR numPhiSteps = root:Packages:NIST:gNPhiSteps713 nphi = numPhiSteps //number of anular sectors is set by users714 715 // rc = 2*dtdist*asin(qc*lambda/4/Pi) //in mm716 // delr = nw*sx/2717 // rlo = rcdelr718 // rhi = rc + delr719 // dphi = 360/nphi720 //721 // /// data wave is data in the current folder which was set at the top of the function722 // Wave data=$(destPath + ":data")723 // //Check for the existence of the mask, if not, make one (local to this folder) that is null724 //725 // if(WaveExists($"root:Packages:NIST:MSK:data") == 0)726 // Print "There is no mask file loaded (WaveExists) the data is not masked"727 // Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")728 // WAVE mask = $(destPath + ":mask")729 // mask = 0730 // else731 // Wave mask=$"root:Packages:NIST:MSK:data"732 // Endif733 //734 // rcentr = 150 //pixels within rcentr of beam center are broken into 9 parts735 // // values for error if unable to estimate value736 // //large_num = 1e10737 // large_num = 1 //1e10 value (typically sig of last data point) plots poorly, arb set to 1738 // small_num = 1e10739 //740 // // output wave are expected to exist (?) initialized to zero, what length?741 // // 300 points on VAX 742 // Variable wavePts=500743 // Make/O/N=(wavePts) $(destPath + ":phival"),$(destPath + ":aveint")744 // Make/O/N=(wavePts) $(destPath + ":ncells"),$(destPath + ":sig"),$(destPath + ":sigave")745 // WAVE phival = $(destPath + ":phival")746 // WAVE aveint = $(destPath + ":aveint")747 // WAVE ncells = $(destPath + ":ncells")748 // WAVE sig = $(destPath + ":sig")749 // WAVE sigave = $(destPath + ":sigave")750 //751 // phival = 0752 // aveint = 0753 // ncells = 0754 // sig = 0755 // sigave = 0756 //757 // dtdis2 = dtdist^2758 // nq = 1759 // xoffst=0760 // //distance of beam center from detector center761 // xbm = FX(x0,sx3,xcenter,sx)762 // ybm = FY(y0,sy3,ycenter,sy)763 //764 // //BEGIN AVERAGE **********765 // Variable xi,xd,x,y,yd,yj,nd,fd,nd2,iphi,ntotal,var766 // Variable jj,data_pixel,xx,yy,ll,kk,rij,phiij,avesq,aveisq767 //768 // // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)769 // // loop index corresponds to FORTRAN (old code)770 // // and the IGOR array indices must be adjusted (1) to the correct address771 // ntotal = 0772 // ii=1773 // do774 // xi = ii775 // xd = FX(xi,sx3,xcenter,sx)776 // x = xoffst + xd xbm //x and y are in mm777 //778 // jj = 1779 // do780 // data_pixel = data[ii1][jj1] //assign to local variable781 // yj = jj782 // yd = FY(yj,sy3,ycenter,sy)783 // y = yd  ybm784 // if(!(mask[ii1][jj1])) //masked pixels = 1, skip if masked (this way works...)785 // nd = 1786 // fd = 1787 // if( (abs(x) > rcentr)  (abs(y) > rcentr)) //break pixel into 9 equal parts788 // nd = 3789 // fd = 2790 // Endif791 // nd2 = nd^2792 // ll = 1 //"elel" loop index793 // do794 // xx = x + (ll  fd)*sx/3795 // kk = 1796 // do797 // yy = y + (kk  fd)*sy/3798 // //test to see if center of pixel (i,j) lies in annulus799 // rij = sqrt(x*x + y*y)/dr + 1.001800 // //check whether pixel lies within width band801 // if((rij > rlo) && (rij < rhi))802 // //in the annulus, do something803 // if (yy >= 0)804 // //phiij is in degrees805 // phiij = atan2(yy,xx)*180/Pi //0 to 180 deg806 // else807 // phiij = 360 + atan2(yy,xx)*180/Pi //180 to 360 deg808 // Endif809 // if (phiij > (3600.5*dphi))810 // phiij = 360811 // Endif812 // iphi = trunc(phiij/dphi + 1.501)813 // aveint[iphi1] += 9*data_pixel/nd2814 // sig[iphi1] += 9*data_pixel*data_pixel/nd2815 // ncells[iphi1] += 9/nd2816 // ntotal += 9/nd2817 // Endif //check if in annulus818 // kk+=1819 // while(kk<=nd)820 // ll += 1821 // while(ll<=nd)822 // Endif //masked pixel check823 // jj += 1824 // while (jj<=pixelsY)825 // ii += 1826 // while(ii<=pixelsX) //end of the averaging827 //828 // //compute phivalues and errors829 //830 // ntotal /=9831 //832 // kk = 1833 // do834 // phival[kk1] = dphi*(kk1)835 // if(ncells[kk1] != 0)836 // aveint[kk1] = aveint[kk1]/ncells[kk1]837 // avesq = aveint[kk1]*aveint[kk1]838 // aveisq = sig[kk1]/ncells[kk1]839 // var = aveisq  avesq840 // if (var <=0 )841 // sig[kk1] = 0842 // sigave[kk1] = 0843 // ncells[kk1] /=9844 // else845 // if(ncells[kk1] > 9)846 // sigave[kk1] = sqrt(9*var/(ncells[kk1]9))847 // sig[kk1] = sqrt( abs(aveint[kk1])/(ncells[kk1]/9) )848 // ncells[kk1] /=9849 // else850 // sig[kk1] = 0851 // sigave[kk1] = 0852 // ncells[kk1] /=9853 // Endif854 // Endif855 // Endif856 // kk+=1857 // while(kk<=nphi)858 //859 // // data waves were defined as 200 points (=wavePts), but now have less than that (nphi) points860 // // use DeletePoints to remove junk from end of waves861 // Variable startElement,numElements862 // startElement = nphi863 // numElements = wavePts  startElement864 // DeletePoints startElement,numElements, phival,aveint,ncells,sig,sigave865 //866 // //////////////end of VAX Phibin.for867 //868 // //angle dependent transmission correction is not done in phiave869 // Ann_1D_Graph(aveint,phival,sigave)870 871 //get rid of the default mask, if one was created (it is in the current folder)872 //don't just kill "mask" since it might be pointing to the one in the MSK folder873 Killwaves/z $(destPath+":mask")874 875 //return to root folder (redundant)876 SetDataFolder root:877 878 Return 0879 End
Note: See TracChangeset
for help on using the changeset viewer.