Ignore:
Timestamp:
Jun 1, 2011 11:30:33 AM (12 years ago)
Author:
srkline
Message:

Smear_2D procedures to plot resolution function at each pixel have been improved, and to integrate the resolution to check that it normalizes to 1.

+ other minor tweaks to convert VAX time and to clean up list behavior.

Location:
sans/Dev/trunk/NCNR_User_Procedures/Common
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtilsMacro_v40.ipf

    r800 r801  
    11871187 
    11881188        // from 2010 model functions 
    1189         list = RemoveFromList("fTwoYukawa;DoBoxGraph;",list,";") 
     1189        list = RemoveFromList("fTwoYukawa;DoBoxGraph;Gauss2D_theta;",list,";") 
    11901190         
    11911191        // from Debye Sphere method 
  • sans/Dev/trunk/NCNR_User_Procedures/Common/Smear_2D.ipf

    r794 r801  
    377377//  
    378378 
    379 //Function PlotResolution_atPixel(type,xx,yy) 
    380 //      String type 
    381 //      Variable xx,yy 
    382 //       
    383 /////from QxQyExport 
    384 //      String destStr="",typeStr="" 
    385 //      Variable step=1,refnum 
    386 //      destStr = "root:Packages:NIST:"+type 
    387 //       
    388 //      //must select the linear_data to export 
    389 //      NVAR isLog = $(destStr+":gIsLogScale") 
    390 //      if(isLog==1) 
    391 //              typeStr = ":linear_data" 
    392 //      else 
    393 //              typeStr = ":data" 
    394 //      endif 
    395 //       
    396 //      NVAR pixelsX = root:myGlobals:gNPixelsX 
    397 //      NVAR pixelsY = root:myGlobals:gNPixelsY 
    398 //       
    399 //      Wave data=$(destStr+typeStr) 
    400 //      WAVE intw=$(destStr + ":integersRead") 
    401 //      WAVE rw=$(destStr + ":realsRead") 
    402 //      WAVE/T textw=$(destStr + ":textRead") 
    403 //       
    404 ////    Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val,phi,r_dist 
    405 //      Variable qx_val,qy_val,z_val,qval,qz_val,phi,r_dist 
    406 // 
    407 //      Variable xctr,yctr,sdd,lambda,pixSize 
    408 //      xctr = rw[16] 
    409 //      yctr = rw[17] 
    410 //      sdd = rw[18] 
    411 //      lambda = rw[26] 
    412 //      pixSize = rw[13]/10             //convert mm to cm (x and y are the same size pixels) 
    413 //       
    414 ////    qx_val = CalcQx(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)          //+1 converts to detector coordinate system 
    415 ////    qy_val = CalcQy(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
    416 //       
    417 //      qx_val = CalcQx(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)                //+1 converts to detector coordinate system 
    418 //      qy_val = CalcQy(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
    419 //       
    420 ////    Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val 
    421 // 
    422 /////************ 
    423 //// do everything to write out the resolution too 
    424 //      // un-comment these if you want to write out qz_val and qval too, then use the proper save command 
    425 //      qval = CalcQval(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
    426 //      qz_val = CalcQz(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
    427 //      phi = FindPhi( pixSize*((xx+1)-xctr) , pixSize*((yy+1)-yctr))           //(dx,dy) 
    428 //      r_dist = sqrt(  (pixSize*((xx+1)-xctr))^2 +  (pixSize*((yy+1)-yctr))^2 )                //radial distance from ctr to pt 
    429 ////    Redimension/N=(pixelsX*pixelsY) qz_val,qval,phi,r_dist 
    430 //      //everything in 1D now 
    431 ////    Duplicate/O qval SigmaQX,SigmaQY,fsubS 
    432 //      Variable SigmaQX,SigmaQY,fsubS 
    433 // 
    434 //      Variable L2 = rw[18] 
    435 //      Variable BS = rw[21] 
    436 //      Variable S1 = rw[23] 
    437 //      Variable S2 = rw[24] 
    438 //      Variable L1 = rw[25] 
    439 //      Variable lambdaWidth = rw[27]    
    440 //      Variable usingLenses = rw[28]           //new 2007 
    441 // 
    442 //      //Two parameters DDET and APOFF are instrument dependent.  Determine 
    443 //      //these from the instrument name in the header. 
    444 //      //From conversation with JB on 01.06.99 these are the current good values 
    445 //      Variable DDet 
    446 //      NVAR apOff = root:myGlobals:apOff               //in cm 
    447 //      DDet = rw[10]/10                        // header value (X) is in mm, want cm here 
    448 // 
    449 //      Variable ret1,ret2,ret3,del_r 
    450 //      del_r = rw[10] 
    451 //       
    452 //      get2DResolution(qval,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,ret1,ret2,ret3) 
    453 //      SigmaQX = ret1   
    454 //      SigmaQY = ret2   
    455 //      fsubs = ret3     
    456 /////// 
    457 // 
    458 ////    Variable theta,phi,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0 
    459 //      Variable theta,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0,maxSig,nStdDev=3,normFactor 
    460 //      Variable qx_ret,qy_ret 
    461 //       
    462 ////    theta = FindPhi(qx_val,qy_val) 
    463 //// need to rotate properly - theta is defined a starting from +y axis, moving CW 
    464 //// we define phi starting from +x and moving CCW 
    465 //      theta = -phi                    //seems to give the right behavior... 
    466 //       
    467 //       
    468 //      Print qx_val,qy_val,qval 
    469 //      Print "phi, theta",phi,theta 
    470 //       
    471 //      FindQxQy(qval,phi,qx_ret,qy_ret) 
    472 //       
    473 //      sx = SigmaQx 
    474 //      sy = sigmaQy 
    475 //      x0 = qx_val 
    476 //      y0 = qy_val 
    477 //       
    478 //      a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy) 
    479 //      b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy) 
    480 //      c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy) 
    481 //       
    482 //      normFactor = pi/sqrt(a*c-b*b) 
    483 // 
    484 //      Make/O/D/N=(num,num) res 
    485 //      // so the resolution function 'looks right' on a 2D plot - otherwise it always looks like a circle 
    486 //      maxSig = max(sx,sy) 
    487 //      Setscale/I x -nStdDev*maxSig+x0,nStdDev*maxSig+x0,res 
    488 //      Setscale/I y -nStdDev*maxSig+y0,nStdDev*maxSig+y0,res 
    489 ////    Setscale/I x -nStdDev*sx+x0,nStdDev*sx+x0,res 
    490 ////    Setscale/I y -nStdDev*sy+y0,nStdDev*sy+y0,res 
    491 //       
    492 //      Variable xPt,yPt,delx,dely,offx,offy 
    493 //      delx = DimDelta(res,0) 
    494 //      dely = DimDelta(res,1) 
    495 //      offx = DimOffset(res,0) 
    496 //      offy = DimOffset(res,1) 
    497 // 
    498 //      Print "sx,sy = ",sx,sy 
    499 //      for(ii=0;ii<num;ii+=1) 
    500 //              xPt = offx + ii*delx 
    501 //              for(jj=0;jj<num;jj+=1) 
    502 //                      yPt = offy + jj*dely 
    503 //                      res[ii][jj] = exp(-1*(a*(xPt-x0)^2 + 2*b*(xPt-x0)*(yPt-y0) + c*(yPt-y0)^2)) 
    504 //              endfor 
    505 //      endfor   
    506 //      res /= normFactor 
    507 //       
    508 //      //Print sum(res,-inf,inf)*delx*dely 
    509 //      if(WaveExists($"coef")==0) 
    510 //              Make/O/D/N=6 coef 
    511 //      endif 
    512 //      Wave coef=coef 
    513 //      coef[0] = 1 
    514 //      coef[1] = qx_val 
    515 //      coef[2] = qy_val 
    516 //      coef[3] = sx 
    517 //      coef[4] = sy 
    518 //      coef[5] = theta 
    519 // 
    520 ////    Variable t1=StopMSTimer(-2) 
    521 // 
    522 //// 
    523 ////    do2dIntegrationGauss(-nStdDev*maxSig+x0,nStdDev*maxSig+x0,-nStdDev*maxSig+y0,nStdDev*maxSig+y0) 
    524 //// 
    525 // 
    526 ////    Variable elap = (StopMSTimer(-2) - t1)/1e6 
    527 ////    Print "elapsed time = ",elap 
    528 ////    Print "time for 16384 = (minutes)",16384*elap/60 
    529 //      return(0) 
    530 //End 
    531  
    532  
    533 //// this is called each time to integrate the gaussian 
    534 //Function do2dIntegrationGauss(xMin,xMax,yMin,yMax) 
    535 //      Variable xMin,xMax,yMin,yMax 
    536 //       
    537 //      Variable/G globalXmin=xMin 
    538 //      Variable/G globalXmax=xMax 
    539 //      Variable/G globalY 
    540 //                       
    541 //      Variable result=Integrate1d(Gauss2DFuncOuter,yMin,yMax,2,5)         
    542 //      KillVariables/z globalXmax,globalXmin,globalY 
    543 //      print "integration of 2D = ",result 
    544 //End 
    545 // 
    546 //Function Gauss2DFuncOuter(inY) 
    547 //      Variable inY 
    548 //       
    549 //      NVAR globalXmin,globalXmax,globalY 
    550 //      globalY=inY 
    551 //       
    552 //      return integrate1D(Gauss2DFuncInner,globalXmin,globalXmax,2,5)           
    553 //End 
    554 // 
    555 //Function Gauss2DFuncInner(inX) 
    556 //      Variable inX 
    557 //       
    558 //      NVAR globalY 
    559 //      Wave coef=coef 
    560 //       
    561 //      return Gauss2D_theta(coef,inX,GlobalY) 
    562 //End 
    563 // 
    564 //Function Gauss2D_theta(w,x,y) 
    565 //      Wave w 
    566 //      Variable x,y 
    567 //       
    568 //      Variable val,a,b,c 
    569 //      Variable scale,x0,y0,sx,sy,theta,normFactor 
    570 //       
    571 //      scale = w[0] 
    572 //      x0 = w[1] 
    573 //      y0 = w[2] 
    574 //      sx = w[3] 
    575 //      sy = w[4] 
    576 //      theta = w[5] 
    577 //       
    578 //      a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy) 
    579 //      b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy) 
    580 //      c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy) 
    581 //       
    582 //      val = exp(-1*(a*(x-x0)^2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)^2)) 
    583 //       
    584 //      normFactor = pi/sqrt(a*c-b*b) 
    585 //       
    586 //      return(scale*val/normFactor) 
    587 //end 
    588 // 
     379Function PlotResolution_atPixel(type,xx,yy) 
     380        String type 
     381        Variable xx,yy 
     382         
     383///from QxQyExport 
     384        String destStr="",typeStr="" 
     385        Variable step=1,refnum 
     386        destStr = "root:Packages:NIST:"+type 
     387         
     388        //must select the linear_data to export 
     389        NVAR isLog = $(destStr+":gIsLogScale") 
     390        if(isLog==1) 
     391                typeStr = ":linear_data" 
     392        else 
     393                typeStr = ":data" 
     394        endif 
     395         
     396        NVAR pixelsX = root:myGlobals:gNPixelsX 
     397        NVAR pixelsY = root:myGlobals:gNPixelsY 
     398         
     399        Wave data=$(destStr+typeStr) 
     400        WAVE intw=$(destStr + ":integersRead") 
     401        WAVE rw=$(destStr + ":realsRead") 
     402        WAVE/T textw=$(destStr + ":textRead") 
     403         
     404//      Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val,phi,r_dist 
     405        Variable qx_val,qy_val,z_val,qval,qz_val,phi,r_dist 
     406 
     407        Variable xctr,yctr,sdd,lambda,pixSize 
     408        xctr = rw[16] 
     409        yctr = rw[17] 
     410        sdd = rw[18] 
     411        lambda = rw[26] 
     412        pixSize = rw[13]/10             //convert mm to cm (x and y are the same size pixels) 
     413         
     414//      qx_val = CalcQx(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)          //+1 converts to detector coordinate system 
     415//      qy_val = CalcQy(p+1,q+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     416         
     417        qx_val = CalcQx(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10)                //+1 converts to detector coordinate system 
     418        qy_val = CalcQy(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     419         
     420//      Redimension/N=(pixelsX*pixelsY) qx_val,qy_val,z_val 
     421 
     422///************ 
     423// do everything to write out the resolution too 
     424        // un-comment these if you want to write out qz_val and qval too, then use the proper save command 
     425        qval = CalcQval(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     426        qz_val = CalcQz(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) 
     427        phi = FindPhi( pixSize*((xx+1)-xctr) , pixSize*((yy+1)-yctr))           //(dx,dy) 
     428        r_dist = sqrt(  (pixSize*((xx+1)-xctr))^2 +  (pixSize*((yy+1)-yctr))^2 )                //radial distance from ctr to pt 
     429//      Redimension/N=(pixelsX*pixelsY) qz_val,qval,phi,r_dist 
     430        //everything in 1D now 
     431//      Duplicate/O qval SigmaQX,SigmaQY,fsubS 
     432        Variable SigmaQX,SigmaQY,fsubS 
     433 
     434        Variable L2 = rw[18] 
     435        Variable BS = rw[21] 
     436        Variable S1 = rw[23] 
     437        Variable S2 = rw[24] 
     438        Variable L1 = rw[25] 
     439        Variable lambdaWidth = rw[27]    
     440        Variable usingLenses = rw[28]           //new 2007 
     441 
     442        //Two parameters DDET and APOFF are instrument dependent.  Determine 
     443        //these from the instrument name in the header. 
     444        //From conversation with JB on 01.06.99 these are the current good values 
     445        Variable DDet 
     446        NVAR apOff = root:myGlobals:apOff               //in cm 
     447        DDet = rw[10]/10                        // header value (X) is in mm, want cm here 
     448 
     449        Variable ret1,ret2,ret3,del_r 
     450        del_r = rw[10] 
     451         
     452        get2DResolution(qval,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,ret1,ret2,ret3) 
     453        SigmaQX = ret1   
     454        SigmaQY = ret2   
     455        fsubs = ret3     
     456///// 
     457 
     458//      Variable theta,phi,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0 
     459        Variable theta,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0,maxSig,nStdDev=3,normFactor 
     460        Variable qx_ret,qy_ret 
     461         
     462//      theta = FindPhi(qx_val,qy_val) 
     463// need to rotate properly - theta is defined a starting from +y axis, moving CW 
     464// we define phi starting from +x and moving CCW 
     465        theta = -phi                    //seems to give the right behavior... 
     466         
     467         
     468        Print qx_val,qy_val,qval 
     469        Print "phi, theta",phi,theta 
     470         
     471        FindQxQy(qval,phi,qx_ret,qy_ret) 
     472         
     473        sx = SigmaQx 
     474        sy = sigmaQy 
     475        x0 = qx_val 
     476        y0 = qy_val 
     477         
     478        a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy) 
     479        b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy) 
     480        c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy) 
     481         
     482        normFactor = pi/sqrt(a*c-b*b) 
     483 
     484        Make/O/D/N=(num,num) res 
     485        // so the resolution function 'looks right' on a 2D plot - otherwise it always looks like a circle 
     486        maxSig = max(sx,sy) 
     487        Setscale/I x -nStdDev*maxSig+x0,nStdDev*maxSig+x0,res 
     488        Setscale/I y -nStdDev*maxSig+y0,nStdDev*maxSig+y0,res 
     489//      Setscale/I x -nStdDev*sx+x0,nStdDev*sx+x0,res 
     490//      Setscale/I y -nStdDev*sy+y0,nStdDev*sy+y0,res 
     491         
     492        Variable xPt,yPt,delx,dely,offx,offy 
     493        delx = DimDelta(res,0) 
     494        dely = DimDelta(res,1) 
     495        offx = DimOffset(res,0) 
     496        offy = DimOffset(res,1) 
     497 
     498        Print "sx,sy = ",sx,sy 
     499        for(ii=0;ii<num;ii+=1) 
     500                xPt = offx + ii*delx 
     501                for(jj=0;jj<num;jj+=1) 
     502                        yPt = offy + jj*dely 
     503                        res[ii][jj] = exp(-1*(a*(xPt-x0)^2 + 2*b*(xPt-x0)*(yPt-y0) + c*(yPt-y0)^2)) 
     504                endfor 
     505        endfor   
     506        res /= normFactor 
     507         
     508        //Print sum(res,-inf,inf)*delx*dely 
     509        if(WaveExists($"coef")==0) 
     510                Make/O/D/N=6 coef 
     511        endif 
     512        Wave coef=coef 
     513        coef[0] = 1 
     514        coef[1] = qx_val 
     515        coef[2] = qy_val 
     516        coef[3] = sx 
     517        coef[4] = sy 
     518        coef[5] = theta 
     519 
     520//      Variable t1=StopMSTimer(-2) 
     521 
     522// 
     523        do2dIntegrationGauss(-nStdDev*maxSig+x0,nStdDev*maxSig+x0,-nStdDev*maxSig+y0,nStdDev*maxSig+y0) 
     524// 
     525 
     526//      Variable elap = (StopMSTimer(-2) - t1)/1e6 
     527//      Print "elapsed time = ",elap 
     528//      Print "time for 16384 = (minutes)",16384*elap/60 
     529        return(0) 
     530End 
     531 
     532 
     533// this is called each time to integrate the gaussian 
     534Function do2dIntegrationGauss(xMin,xMax,yMin,yMax) 
     535        Variable xMin,xMax,yMin,yMax 
     536         
     537        Variable/G globalXmin=xMin 
     538        Variable/G globalXmax=xMax 
     539        Variable/G globalY 
     540                         
     541        Variable result=Integrate1d(Gauss2DFuncOuter,yMin,yMax,2,5)         
     542        KillVariables/z globalXmax,globalXmin,globalY 
     543        print "integration of 2D = ",result 
     544End 
     545 
     546Function Gauss2DFuncOuter(inY) 
     547        Variable inY 
     548         
     549        NVAR globalXmin,globalXmax,globalY 
     550        globalY=inY 
     551         
     552        return integrate1D(Gauss2DFuncInner,globalXmin,globalXmax,2,5)           
     553End 
     554 
     555Function Gauss2DFuncInner(inX) 
     556        Variable inX 
     557         
     558        NVAR globalY 
     559        Wave coef=coef 
     560         
     561        return Gauss2D_theta(coef,inX,GlobalY) 
     562End 
     563 
     564Function Gauss2D_theta(w,x,y) 
     565        Wave w 
     566        Variable x,y 
     567         
     568        Variable val,a,b,c 
     569        Variable scale,x0,y0,sx,sy,theta,normFactor 
     570         
     571        scale = w[0] 
     572        x0 = w[1] 
     573        y0 = w[2] 
     574        sx = w[3] 
     575        sy = w[4] 
     576        theta = w[5] 
     577         
     578        a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy) 
     579        b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy) 
     580        c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy) 
     581         
     582        val = exp(-1*(a*(x-x0)^2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)^2)) 
     583         
     584        normFactor = pi/sqrt(a*c-b*b) 
     585         
     586        return(scale*val/normFactor) 
     587end 
     588 
    589589//////////////////////////////////////////////////////////// 
    590590//////////////////////////////////////////////////////////// 
Note: See TracChangeset for help on using the changeset viewer.