Changeset 801
- Timestamp:
- Jun 1, 2011 11:30:33 AM (12 years ago)
- Location:
- sans/Dev/trunk/NCNR_User_Procedures
- Files:
-
- 5 edited
-
Analysis/Packages/Wrapper_v40.ipf (modified) (2 diffs)
-
Common/Packages/PlotManager/PlotUtilsMacro_v40.ipf (modified) (1 diff)
-
Common/Smear_2D.ipf (modified) (1 diff)
-
Reduction/SANS/Initialize.ipf (modified) (1 diff)
-
Reduction/SANS/NCNR_Utils.ipf (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
sans/Dev/trunk/NCNR_User_Procedures/Analysis/Packages/Wrapper_v40.ipf
r764 r801 942 942 // Variable t0 = stopMStimer(-2) 943 943 944 FuncFit/H=getHStr(hold) / NTHR=0 $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /STRC=fs /NWOK944 FuncFit/H=getHStr(hold) /M=2 /NTHR=0 $funcStr cw, yw[pt1,pt2] /X=xw /W=sw /I=1 /E=eps /D=fitYw /C=constr /STRC=fs /NWOK 945 945 946 946 // t0 = (stopMSTimer(-2) - t0)*1e-6 … … 1054 1054 1055 1055 if(yesReport) 1056 String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+ suffix, "", "TEXT:1," ) // this is *hopefully* one wave1056 String parStr=GetWavesDataFolder(cw,1)+ WaveList("*param*"+"_"+suffix, "", "TEXT:1," ) // this is *hopefully* one wave 1057 1057 String topGraph= WinName(0,1) //this is the topmost graph 1058 1058 -
sans/Dev/trunk/NCNR_User_Procedures/Common/Packages/PlotManager/PlotUtilsMacro_v40.ipf
r800 r801 1187 1187 1188 1188 // from 2010 model functions 1189 list = RemoveFromList("fTwoYukawa;DoBoxGraph; ",list,";")1189 list = RemoveFromList("fTwoYukawa;DoBoxGraph;Gauss2D_theta;",list,";") 1190 1190 1191 1191 // from Debye Sphere method -
sans/Dev/trunk/NCNR_User_Procedures/Common/Smear_2D.ipf
r794 r801 377 377 // 378 378 379 //Function PlotResolution_atPixel(type,xx,yy)380 //String type381 //Variable xx,yy382 //383 /// //from QxQyExport384 //String destStr="",typeStr=""385 //Variable step=1,refnum386 //destStr = "root:Packages:NIST:"+type387 //388 ////must select the linear_data to export389 //NVAR isLog = $(destStr+":gIsLogScale")390 //if(isLog==1)391 //typeStr = ":linear_data"392 //else393 //typeStr = ":data"394 //endif395 //396 //NVAR pixelsX = root:myGlobals:gNPixelsX397 //NVAR pixelsY = root:myGlobals:gNPixelsY398 //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_dist405 //Variable qx_val,qy_val,z_val,qval,qz_val,phi,r_dist406 // 407 //Variable xctr,yctr,sdd,lambda,pixSize408 //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 system415 // //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 system418 //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_val421 // 422 /// //************423 // //do everything to write out the resolution too424 //// un-comment these if you want to write out qz_val and qval too, then use the proper save command425 //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 pt429 // //Redimension/N=(pixelsX*pixelsY) qz_val,qval,phi,r_dist430 ////everything in 1D now431 // //Duplicate/O qval SigmaQX,SigmaQY,fsubS432 //Variable SigmaQX,SigmaQY,fsubS433 // 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 2007441 // 442 ////Two parameters DDET and APOFF are instrument dependent. Determine443 ////these from the instrument name in the header.444 ////From conversation with JB on 01.06.99 these are the current good values445 //Variable DDet446 //NVAR apOff = root:myGlobals:apOff //in cm447 //DDet = rw[10]/10 // header value (X) is in mm, want cm here448 // 449 //Variable ret1,ret2,ret3,del_r450 //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 = ret1454 //SigmaQY = ret2455 //fsubs = ret3456 ///// //457 // 458 // //Variable theta,phi,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0459 //Variable theta,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0,maxSig,nStdDev=3,normFactor460 //Variable qx_ret,qy_ret461 //462 // //theta = FindPhi(qx_val,qy_val)463 // //need to rotate properly - theta is defined a starting from +y axis, moving CW464 // //we define phi starting from +x and moving CCW465 //theta = -phi //seems to give the right behavior...466 //467 //468 //Print qx_val,qy_val,qval469 //Print "phi, theta",phi,theta470 //471 //FindQxQy(qval,phi,qx_ret,qy_ret)472 //473 //sx = SigmaQx474 //sy = sigmaQy475 //x0 = qx_val476 //y0 = qy_val477 //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) res485 //// so the resolution function 'looks right' on a 2D plot - otherwise it always looks like a circle486 //maxSig = max(sx,sy)487 //Setscale/I x -nStdDev*maxSig+x0,nStdDev*maxSig+x0,res488 //Setscale/I y -nStdDev*maxSig+y0,nStdDev*maxSig+y0,res489 // //Setscale/I x -nStdDev*sx+x0,nStdDev*sx+x0,res490 // //Setscale/I y -nStdDev*sy+y0,nStdDev*sy+y0,res491 //492 //Variable xPt,yPt,delx,dely,offx,offy493 //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,sy499 //for(ii=0;ii<num;ii+=1)500 //xPt = offx + ii*delx501 //for(jj=0;jj<num;jj+=1)502 //yPt = offy + jj*dely503 //res[ii][jj] = exp(-1*(a*(xPt-x0)^2 + 2*b*(xPt-x0)*(yPt-y0) + c*(yPt-y0)^2))504 //endfor505 //endfor506 //res /= normFactor507 //508 ////Print sum(res,-inf,inf)*delx*dely509 //if(WaveExists($"coef")==0)510 //Make/O/D/N=6 coef511 //endif512 //Wave coef=coef513 //coef[0] = 1514 //coef[1] = qx_val515 //coef[2] = qy_val516 //coef[3] = sx517 //coef[4] = sy518 //coef[5] = theta519 // 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)/1e6527 // //Print "elapsed time = ",elap528 // //Print "time for 16384 = (minutes)",16384*elap/60529 //return(0)530 //End531 532 533 // //this is called each time to integrate the gaussian534 //Function do2dIntegrationGauss(xMin,xMax,yMin,yMax)535 //Variable xMin,xMax,yMin,yMax536 //537 //Variable/G globalXmin=xMin538 //Variable/G globalXmax=xMax539 //Variable/G globalY540 //541 //Variable result=Integrate1d(Gauss2DFuncOuter,yMin,yMax,2,5)542 //KillVariables/z globalXmax,globalXmin,globalY543 //print "integration of 2D = ",result544 //End545 // 546 //Function Gauss2DFuncOuter(inY)547 //Variable inY548 //549 //NVAR globalXmin,globalXmax,globalY550 //globalY=inY551 //552 //return integrate1D(Gauss2DFuncInner,globalXmin,globalXmax,2,5)553 //End554 // 555 //Function Gauss2DFuncInner(inX)556 //Variable inX557 //558 //NVAR globalY559 //Wave coef=coef560 //561 //return Gauss2D_theta(coef,inX,GlobalY)562 //End563 // 564 //Function Gauss2D_theta(w,x,y)565 //Wave w566 //Variable x,y567 //568 //Variable val,a,b,c569 //Variable scale,x0,y0,sx,sy,theta,normFactor570 //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 //end588 // 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 589 589 //////////////////////////////////////////////////////////// 590 590 //////////////////////////////////////////////////////////// -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/Initialize.ipf
r749 r801 63 63 Main_Panel() 64 64 Endif 65 ResizeCmdWindow()65 // ResizeCmdWindow() 66 66 67 67 //unload the NCNR_Package_Loader, if NCNR not defined -
sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/NCNR_Utils.ipf
r800 r801 520 520 // print yr,mon,day 521 521 time_secs = date2secs(yr,mon,day) 522 523 return(time_secs) 524 end 525 526 // dd-mon-yyyy hh:mm:ss -> seconds 527 // the VAX uses 24 hr time for hh 528 // 529 Function ConvertVAXDateTime2secs(dateAndTime) 530 string dateAndTime 531 532 Variable day,yr,mon,hh,mm,ss,time_secs 533 string str,monStr 534 535 str=dateandtime 536 sscanf str,"%d-%3s-%4d %d:%d:%d",day,monStr,yr,hh,mm,ss 537 mon = monStr2num(monStr) 538 // print yr,mon,day,hh,mm,ss 539 time_secs = date2secs(yr,mon,day) 540 time_secs += hh*3600 + mm*60 + ss 541 522 542 523 543 return(time_secs) … … 1501 1521 End 1502 1522 1523 // a utility function so that I can get the values from the command line 1524 // since the atten_err is PBR 1525 // 1526 Function PrintAttenuation(instr,lam,attenNo) 1527 String instr 1528 Variable lam,attenNo 1529 1530 Variable atten_err, attenFactor 1531 1532 strswitch(instr) 1533 case "NG3": 1534 attenFactor = LookupAttenNG3(lam,attenNo,atten_err) 1535 break 1536 case "NG5": 1537 //using NG7 lookup Table 1538 attenFactor = LookupAttenNG7(lam,attenNo,atten_err) 1539 break 1540 case "NG7": 1541 attenFactor = LookupAttenNG7(lam,attenNo,atten_err) 1542 break 1543 default: 1544 //return error? 1545 attenFactor=1 1546 endswitch 1547 1548 Print "atten, err = ", attenFactor, atten_err 1549 1550 return(0) 1551 End 1552 1553 1554 1503 1555 //returns the proper attenuation factor based on the instrument (NG3, NG5, or NG7) 1504 1556 //NG5 values are taken from the NG7 tables (there is very little difference in the
Note: See TracChangeset
for help on using the changeset viewer.
