Changeset 1208


Ignore:
Timestamp:
Sep 13, 2019 3:19:52 PM (3 years ago)
Author:
krzywon
Message:

Add an ellipse drawing tool that (mostly) works, and start on the elliptical averaging algorithm.

Location:
sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/AvgGraphics.ipf

    r1207 r1208  
    327327        GroupBox ann,pos={148,44},size={143,84},title="Annular" 
    328328        GroupBox rect,pos={7,133},size={134,71},title="Rectangular" 
    329         GroupBox sect,pos={148,133},size={143,71},title="Sector" 
     329        GroupBox sect,pos={148,133},size={143,71},title="Sector/Elliptical" 
    330330        GroupBox sect_rect,pos={7,44},size={134,84},title="Sector/Rectangular" 
    331331        PopupMenu av_choice,pos={61,7},size={144,20},proc=AvTypePopMenuProc,title="AverageType" 
     
    476476                Return 0                //exit the Draw routine 
    477477        Endif 
     478         
     479        // Elliptical 
     480        if(cmpstr(av_type,"Elliptical")==0) 
     481                // 
     482                // TODO: Write the elliptical averaging design 
     483                // 
     484                 
     485                //do the elliptical drawing 
     486                sdd=reals[18] 
     487                lam=reals[26] 
     488                 
     489                sdd *=100               //convert to cm 
     490                //find the radius (from the y-direction) 
     491                variable RAxes = NumberByKey("RATIOAXES",drawStr,"=",";") 
     492                 
     493                Make/O/N=(128,128) ellipsewave 
     494                variable b = 64 
     495                variable a = RAxes * b 
     496                DrawAnEllipsoid("ellipsewave",x0,y0,a,b,phi) 
     497          
     498                Return 0                //exit the Draw routine 
     499        EndIf 
    478500         
    479501        //else sector or rectangular - draw the lines 
     
    795817End 
    796818 
     819  
     820// This creates an ellipsoid 
     821// INPUTS: 
     822//   wname - name of xywave 
     823//   (xo,yo) - center position of ellipse 
     824//   (a, b) - elliptical dimesions to above 
     825//   phi - angle of rotation (counterclockwise) 
     826// OUTPUTS: 
     827//   make all changes directly in xywave 
     828Function DrawAnEllipsoid(wname,xo,yo,a,b,phi) 
     829    string wname 
     830    variable xo,yo,a,b,phi 
     831    int ii 
     832 
     833    wave xywave = $wname 
     834     
     835    make/FREE/n=361 urx,ury 
     836    make/n=361 xw,yw = 0 
     837    
     838    variable cosalpha = cos(phi*Pi/180), sinalpha = sin(phi*Pi/180) 
     839    
     840    urx = a*(cos(p*Pi/180)) 
     841    ury = b*(sin(p*Pi/180)) 
     842    
     843    xw = xo + cosalpha* urx - sinalpha*ury 
     844    yw = yo + sinalpha*urx + cosalpha*ury 
     845     
     846    // FIXME: No single point mapping in Igor - need to do this automatically 
     847     
     848    xywave[xw,yw] = 1 
     849     
     850    rr=0 
     851    gg=50000 
     852    bb=0 
     853     
     854         //SetDrawEnv/W=SANS_Data xcoord= bottom,ycoord= left,linefgc= (rr,gg,bb),linethick= (thick),fillpat=0 
     855         //AppendImage/W=SANS_Data ellipsewave 
     856        
     857    Return (0)  
     858End 
     859 
    797860//reads the value from the panel and updates the global string 
    798861// redraws the appropriate line on the image 
     
    866929         
    867930        String/G root:myGlobals:Drawing:gDrawInfoStr = newStr 
     931         
     932        //redraw the angles 
     933        MasterAngleDraw() 
     934End 
     935 
     936//reads the value from the panel and updates the global string 
     937// redraws the appropriate line on the image 
     938Function AxisRatioSetVarProc(ctrlName,varNum,varStr,varName) : SetVariableControl 
     939        String ctrlName 
     940        Variable varNum 
     941        String varStr 
     942        String varName 
     943 
     944        ControlInfo/W=Average_Panel RAxes_p 
     945        Variable val = V_Value 
     946        SVAR tempStr = root:myGlobals:Drawing:gDrawRAxes 
     947        String newStr = ReplaceNumberByKey("RATIOAXES", tempStr, val, "=", ";") 
     948         
     949        String/G root:myGlobals:Drawing:gDrawRAxes = newStr 
    868950         
    869951        //redraw the angles 
  • sans/Dev/branches/elliptical_averaging/NCNR_User_Procedures/Reduction/SANS/RectAnnulAvg.ipf

    r1207 r1208  
    706706        lambda = reals[26] 
    707707         
    708         Variable qc = NumberByKey("QCENTER",keyListStr,"=",";") 
    709         Variable nw = NumberByKey("QDELTA",keyListStr,"=",";") 
     708        Variable phi = NumberByKey("PHI",keyListStr,"=",";") 
     709        Variable ratio = NumberByKey("RATIOAXES",keyListStr,"=",";") 
    710710         
    711711        dr = 1                  //minimum annulus width, keep this fixed at one 
     
    713713        nphi = numPhiSteps              //number of anular sectors is set by users 
    714714         
    715         rc = 2*dtdist*asin(qc*lambda/4/Pi)              //in mm 
    716         delr = nw*sx/2 
    717         rlo = rc-delr 
    718         rhi = rc + delr 
    719         dphi = 360/nphi 
    720  
    721         /// data wave is data in the current folder which was set at the top of the function 
    722         Wave data=$(destPath + ":data") 
    723         //Check for the existence of the mask, if not, make one (local to this folder) that is null 
    724          
    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 = 0 
    730         else 
    731                 Wave mask=$"root:Packages:NIST:MSK:data" 
    732         Endif 
    733          
    734         rcentr = 150            //pixels within rcentr of beam center are broken into 9 parts 
    735         // values for error if unable to estimate value 
    736         //large_num = 1e10 
    737         large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1 
    738         small_num = 1e-10 
    739          
    740         // output wave are expected to exist (?) initialized to zero, what length? 
    741         // 300 points on VAX --- 
    742         Variable wavePts=500 
    743         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 = 0 
    752         aveint = 0 
    753         ncells = 0 
    754         sig = 0 
    755         sigave = 0 
    756  
    757         dtdis2 = dtdist^2 
    758         nq = 1 
    759         xoffst=0 
    760         //distance of beam center from detector center 
    761         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,var 
    766         Variable jj,data_pixel,xx,yy,ll,kk,rij,phiij,avesq,aveisq 
    767  
    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 address 
    771         ntotal = 0 
    772         ii=1 
    773         do 
    774                 xi = ii 
    775                 xd = FX(xi,sx3,xcenter,sx) 
    776                 x = xoffst + xd -xbm            //x and y are in mm 
    777                  
    778                 jj = 1 
    779                 do 
    780                         data_pixel = data[ii-1][jj-1]           //assign to local variable 
    781                         yj = jj 
    782                         yd = FY(yj,sy3,ycenter,sy) 
    783                         y = yd - ybm 
    784                         if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...) 
    785                                 nd = 1 
    786                                 fd = 1 
    787                                 if( (abs(x) > rcentr) || (abs(y) > rcentr))     //break pixel into 9 equal parts 
    788                                         nd = 3 
    789                                         fd = 2 
    790                                 Endif 
    791                                 nd2 = nd^2 
    792                                 ll = 1          //"el-el" loop index 
    793                                 do 
    794                                         xx = x + (ll - fd)*sx/3 
    795                                         kk = 1 
    796                                         do 
    797                                                 yy = y + (kk - fd)*sy/3 
    798                                                 //test to see if center of pixel (i,j) lies in annulus 
    799                                                 rij = sqrt(x*x + y*y)/dr + 1.001 
    800                                                 //check whether pixel lies within width band 
    801                                                 if((rij > rlo) && (rij < rhi)) 
    802                                                         //in the annulus, do something 
    803                                                         if (yy >= 0) 
    804                                                                 //phiij is in degrees 
    805                                                                 phiij = atan2(yy,xx)*180/Pi             //0 to 180 deg 
    806                                                         else 
    807                                                                 phiij = 360 + atan2(yy,xx)*180/Pi               //180 to 360 deg 
    808                                                         Endif 
    809                                                         if (phiij > (360-0.5*dphi)) 
    810                                                                 phiij -= 360 
    811                                                         Endif 
    812                                                         iphi = trunc(phiij/dphi + 1.501) 
    813                                                         aveint[iphi-1] += 9*data_pixel/nd2 
    814                                                         sig[iphi-1] += 9*data_pixel*data_pixel/nd2 
    815                                                         ncells[iphi-1] += 9/nd2 
    816                                                         ntotal += 9/nd2 
    817                                                 Endif           //check if in annulus 
    818                                                 kk+=1 
    819                                         while(kk<=nd) 
    820                                         ll += 1 
    821                                 while(ll<=nd) 
    822                         Endif           //masked pixel check 
    823                         jj += 1 
    824                 while (jj<=pixelsY) 
    825                 ii += 1 
    826         while(ii<=pixelsX)              //end of the averaging 
    827                  
    828         //compute phi-values and errors 
    829          
    830         ntotal /=9 
    831          
    832         kk = 1 
    833         do 
    834                 phival[kk-1] = dphi*(kk-1) 
    835                 if(ncells[kk-1] != 0) 
    836                         aveint[kk-1] = aveint[kk-1]/ncells[kk-1] 
    837                         avesq = aveint[kk-1]*aveint[kk-1] 
    838                         aveisq = sig[kk-1]/ncells[kk-1] 
    839                         var = aveisq - avesq 
    840                         if (var <=0 ) 
    841                                 sig[kk-1] = 0 
    842                                 sigave[kk-1] = 0 
    843                                 ncells[kk-1] /=9 
    844                         else 
    845                                 if(ncells[kk-1] > 9) 
    846                                         sigave[kk-1] = sqrt(9*var/(ncells[kk-1]-9)) 
    847                                         sig[kk-1] = sqrt( abs(aveint[kk-1])/(ncells[kk-1]/9) ) 
    848                                         ncells[kk-1] /=9 
    849                                 else 
    850                                         sig[kk-1] = 0 
    851                                         sigave[kk-1] = 0 
    852                                         ncells[kk-1] /=9 
    853                                 Endif 
    854                         Endif 
    855                 Endif 
    856                 kk+=1 
    857         while(kk<=nphi) 
    858          
    859         // data waves were defined as 200 points (=wavePts), but now have less than that (nphi) points 
    860         // use DeletePoints to remove junk from end of waves 
    861         Variable startElement,numElements 
    862         startElement = nphi 
    863         numElements = wavePts - startElement 
    864         DeletePoints startElement,numElements, phival,aveint,ncells,sig,sigave 
    865          
    866         //////////////end of VAX Phibin.for 
    867                  
    868         //angle dependent transmission correction is not done in phiave 
    869         Ann_1D_Graph(aveint,phival,sigave) 
     715//      rc = 2*dtdist*asin(qc*lambda/4/Pi)              //in mm 
     716//      delr = nw*sx/2 
     717//      rlo = rc-delr 
     718//      rhi = rc + delr 
     719//      dphi = 360/nphi 
     720// 
     721//      /// data wave is data in the current folder which was set at the top of the function 
     722//      Wave data=$(destPath + ":data") 
     723//      //Check for the existence of the mask, if not, make one (local to this folder) that is null 
     724//       
     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 = 0 
     730//      else 
     731//              Wave mask=$"root:Packages:NIST:MSK:data" 
     732//      Endif 
     733//       
     734//      rcentr = 150            //pixels within rcentr of beam center are broken into 9 parts 
     735//      // values for error if unable to estimate value 
     736//      //large_num = 1e10 
     737//      large_num = 1           //1e10 value (typically sig of last data point) plots poorly, arb set to 1 
     738//      small_num = 1e-10 
     739//       
     740//      // output wave are expected to exist (?) initialized to zero, what length? 
     741//      // 300 points on VAX --- 
     742//      Variable wavePts=500 
     743//      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 = 0 
     752//      aveint = 0 
     753//      ncells = 0 
     754//      sig = 0 
     755//      sigave = 0 
     756// 
     757//      dtdis2 = dtdist^2 
     758//      nq = 1 
     759//      xoffst=0 
     760//      //distance of beam center from detector center 
     761//      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,var 
     766//      Variable jj,data_pixel,xx,yy,ll,kk,rij,phiij,avesq,aveisq 
     767// 
     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 address 
     771//      ntotal = 0 
     772//      ii=1 
     773//      do 
     774//              xi = ii 
     775//              xd = FX(xi,sx3,xcenter,sx) 
     776//              x = xoffst + xd -xbm            //x and y are in mm 
     777//               
     778//              jj = 1 
     779//              do 
     780//                      data_pixel = data[ii-1][jj-1]           //assign to local variable 
     781//                      yj = jj 
     782//                      yd = FY(yj,sy3,ycenter,sy) 
     783//                      y = yd - ybm 
     784//                      if(!(mask[ii-1][jj-1]))                 //masked pixels = 1, skip if masked (this way works...) 
     785//                              nd = 1 
     786//                              fd = 1 
     787//                              if( (abs(x) > rcentr) || (abs(y) > rcentr))     //break pixel into 9 equal parts 
     788//                                      nd = 3 
     789//                                      fd = 2 
     790//                              Endif 
     791//                              nd2 = nd^2 
     792//                              ll = 1          //"el-el" loop index 
     793//                              do 
     794//                                      xx = x + (ll - fd)*sx/3 
     795//                                      kk = 1 
     796//                                      do 
     797//                                              yy = y + (kk - fd)*sy/3 
     798//                                              //test to see if center of pixel (i,j) lies in annulus 
     799//                                              rij = sqrt(x*x + y*y)/dr + 1.001 
     800//                                              //check whether pixel lies within width band 
     801//                                              if((rij > rlo) && (rij < rhi)) 
     802//                                                      //in the annulus, do something 
     803//                                                      if (yy >= 0) 
     804//                                                              //phiij is in degrees 
     805//                                                              phiij = atan2(yy,xx)*180/Pi             //0 to 180 deg 
     806//                                                      else 
     807//                                                              phiij = 360 + atan2(yy,xx)*180/Pi               //180 to 360 deg 
     808//                                                      Endif 
     809//                                                      if (phiij > (360-0.5*dphi)) 
     810//                                                              phiij -= 360 
     811//                                                      Endif 
     812//                                                      iphi = trunc(phiij/dphi + 1.501) 
     813//                                                      aveint[iphi-1] += 9*data_pixel/nd2 
     814//                                                      sig[iphi-1] += 9*data_pixel*data_pixel/nd2 
     815//                                                      ncells[iphi-1] += 9/nd2 
     816//                                                      ntotal += 9/nd2 
     817//                                              Endif           //check if in annulus 
     818//                                              kk+=1 
     819//                                      while(kk<=nd) 
     820//                                      ll += 1 
     821//                              while(ll<=nd) 
     822//                      Endif           //masked pixel check 
     823//                      jj += 1 
     824//              while (jj<=pixelsY) 
     825//              ii += 1 
     826//      while(ii<=pixelsX)              //end of the averaging 
     827//               
     828//      //compute phi-values and errors 
     829//       
     830//      ntotal /=9 
     831//       
     832//      kk = 1 
     833//      do 
     834//              phival[kk-1] = dphi*(kk-1) 
     835//              if(ncells[kk-1] != 0) 
     836//                      aveint[kk-1] = aveint[kk-1]/ncells[kk-1] 
     837//                      avesq = aveint[kk-1]*aveint[kk-1] 
     838//                      aveisq = sig[kk-1]/ncells[kk-1] 
     839//                      var = aveisq - avesq 
     840//                      if (var <=0 ) 
     841//                              sig[kk-1] = 0 
     842//                              sigave[kk-1] = 0 
     843//                              ncells[kk-1] /=9 
     844//                      else 
     845//                              if(ncells[kk-1] > 9) 
     846//                                      sigave[kk-1] = sqrt(9*var/(ncells[kk-1]-9)) 
     847//                                      sig[kk-1] = sqrt( abs(aveint[kk-1])/(ncells[kk-1]/9) ) 
     848//                                      ncells[kk-1] /=9 
     849//                              else 
     850//                                      sig[kk-1] = 0 
     851//                                      sigave[kk-1] = 0 
     852//                                      ncells[kk-1] /=9 
     853//                              Endif 
     854//                      Endif 
     855//              Endif 
     856//              kk+=1 
     857//      while(kk<=nphi) 
     858//       
     859//      // data waves were defined as 200 points (=wavePts), but now have less than that (nphi) points 
     860//      // use DeletePoints to remove junk from end of waves 
     861//      Variable startElement,numElements 
     862//      startElement = nphi 
     863//      numElements = wavePts - startElement 
     864//      DeletePoints startElement,numElements, phival,aveint,ncells,sig,sigave 
     865//       
     866//      //////////////end of VAX Phibin.for 
     867//               
     868//      //angle dependent transmission correction is not done in phiave 
     869//      Ann_1D_Graph(aveint,phival,sigave) 
    870870         
    871871        //get rid of the default mask, if one was created (it is in the current folder) 
Note: See TracChangeset for help on using the changeset viewer.