 Timestamp:
 Feb 3, 2009 3:57:08 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf
r465 r472 56 56 //  a number of models are now found to be illbehaved when q=1e10. Then the random deviate calculation blows up. 57 57 // a warning has been added  but some models need a proper limiting value, and some (powerlaw) are simply unuseable 58 // unless something else can be done. 58 // unless something else can be done. Using a logspacing of points doesn't seem to help, and it introduces a lot of 59 // other problems. Not the way to go. 59 60 //  if the MC gags on a simulation, it often gets "stuck" and can't do the normal calculation from the model, which it 60 61 // should always default to... … … 70 71 // does nothing in the Igor code 71 72 Duplicate/O results retWave 72 //results[0] = 1*(datetime)73 73 74 74 Variable NNeutron=inputWave[0] … … 81 81 82 82 variable mt= ThreadGroupCreate(nthreads) 83 NVAR gInitTime = root:Packages:NIST:SAS:gRanDateTime //time that SASCALC was started 84 83 85 84 86 inputWave[0] = NNeutron/nthreads //split up the number of neutrons … … 93 95 Duplicate/O inputWave $("inputWave"+num2istr(i)) 94 96 Duplicate/O ran_dev $("ran_dev"+num2istr(i)) 95 97 96 98 // ?? I need explicit wave references? 97 99 // maybe I need to have everything in separate data folders  bu tI haven't tried that. seems like a reach. … … 99 101 if(i==0) 100 102 WAVE inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0 101 retWave0[0] = 1*datetime //to initialize ran3 103 retWave0 = 0 //clear the return wave 104 retWave0[0] = 1*(datetimegInitTime) //to initialize ran3 102 105 ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0) 103 106 Print "started thread 0" … … 105 108 if(i==1) 106 109 WAVE inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1 107 //retWave1[0] = 1*datetime //to initialize ran3 110 retWave1 = 0 //clear the return wave 111 retWave1[0] = 1*(datetimegInitTime) //to initialize ran1 108 112 ThreadStart mt,i,Monte_SANS_W2(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1) 109 113 Print "started thread 1" … … 180 184 181 185 // worker function for threads, does nothing except switch between XOP and Igor versions 186 // 187 // uses ran3 182 188 ThreadSafe Function Monte_SANS_W1(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 183 189 WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results … … 191 197 return (0) 192 198 End 199 193 200 // worker function for threads, does nothing except switch between XOP and Igor versions 201 // 202 // uses ran1 194 203 ThreadSafe Function Monte_SANS_W2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 195 204 WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results … … 569 578 qu = 4*pi/lam 570 579 571 // Make/O/N=(nPts_ran)/D root:Packages:NIST:SAS:Gq,root:Packages:NIST:SAS:xw // if these waves are 1000 pts, the results are "pixelated"572 // WAVE Gq = root:Packages:NIST:SAS:gQ573 // WAVE xw = root:Packages:NIST:SAS:xw574 575 580 // hardwired into the Simulation directory rather than the SAS folder. 576 581 // plotting resolutionsmeared models won't work any other way … … 578 583 WAVE Gq = root:Simulation:gQ 579 584 WAVE xw = root:Simulation:xw 580 SetScale/I x (0+1e6),qu*(11e10),"", Gq,xw //don't start at zero or run up all the way to qu to avoid numerical errors 581 582 /// 585 SetScale/I x (0+1e4),qu*(11e10),"", Gq,xw //don't start at zero or run up all the way to qu to avoid numerical errors 586 583 587 /// if all of the coefficients are wellbehaved, then the last point is the background 584 588 // and I can set it to zero here (only for the calculation) … … 593 597 Gq = Gq*sin(2*asin(x/qu))/sqrt(1(x/qu)) // exact 594 598 // 599 // 595 600 Integrate/METH=1 Gq/D=Gq_INT 596 601 … … 605 610 End 606 611 607 612 // returns the random deviate as a wave 613 // and the total SAS crosssection [1/cm] sig_sas 614 // 615 // uses a log spacing of x for better coverage 616 // downside is that it doesn't use builtin integrate, which is automatically cumulative 617 // 618 //  Currently does not work  since the usage of the random deviate in the MC routine is based on the 619 // wave properties of ran_dev, that is it must have the proper scaling and be equally spaced. 620 // 621 //  not really sure that this will solve any of the problems with some functions (notably those with powerlaws) 622 // giving unreasonably large SAS cross sections. (>>10) 623 // 624 Function CalculateRandomDeviate_log(func,coef,lam,outWave,SASxs) 625 FUNCREF SANSModelAAO_MCproto func 626 WAVE coef 627 Variable lam 628 String outWave 629 Variable &SASxs 630 631 Variable nPts_ran=1000,qu,qmin,ii 632 qmin=1e5 633 qu = 4*pi/lam 634 635 // hardwired into the Simulation directory rather than the SAS folder. 636 // plotting resolutionsmeared models won't work any other way 637 Make/O/N=(nPts_ran)/D root:Simulation:Gq,root:Simulation:xw // if these waves are 1000 pts, the results are "pixelated" 638 WAVE Gq = root:Simulation:gQ 639 WAVE xw = root:Simulation:xw 640 // SetScale/I x (0+1e4),qu*(11e10),"", Gq,xw //don't start at zero or run up all the way to qu to avoid numerical errors 641 xw = alog(log(qmin) + x*((log(qu)log(qmin))/nPts_ran)) 642 643 /// if all of the coefficients are wellbehaved, then the last point is the background 644 // and I can set it to zero here (only for the calculation) 645 Duplicate/O coef,tmp_coef 646 Variable num=numpnts(coef) 647 tmp_coef[num1] = 0 648 649 func(tmp_coef,Gq,xw) //call as AAO 650 Gq = Gq*sin(2*asin(xw/qu))/sqrt(1(xw/qu)) // exact 651 652 653 Duplicate/O Gq Gq_INT 654 Gq_INT = 0 655 for(ii=0;ii<nPts_ran;ii+=1) 656 Gq_INT[ii] = AreaXY(xw,Gq,qmin,xw[ii]) 657 endfor 658 659 SASxs = lam*Gq_INT[nPts_ran1] 660 661 Gq_INT /= Gq_INT[nPts_ran1] 662 663 Duplicate/O Gq_INT $outWave 664 665 return(0) 666 End 608 667 609 668 ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) … … 802 861 803 862 // globals are initialized in SASCALC.ipf 863 // coordinates if I make this part of the panel  but this breaks other things... 864 // 865 //Proc MC_SASCALC() 866 //// PauseUpdate; Silent 1 // building window... 867 // 868 //// NewPanel /W=(92,556,390,1028)/K=1 as "SANS Simulator" 869 // SetVariable MC_setvar0,pos={491,73},size={144,15},bodyWidth=80,title="# of neutrons" 870 // SetVariable MC_setvar0,format="%5.4g" 871 // SetVariable MC_setvar0,limits={inf,inf,100},value= root:Packages:NIST:SAS:gImon 872 // SetVariable MC_setvar0_1,pos={491,119},size={131,15},bodyWidth=60,title="Thickness (cm)" 873 // SetVariable MC_setvar0_1,limits={inf,inf,0.1},value= root:Packages:NIST:SAS:gThick 874 // SetVariable MC_setvar0_2,pos={491,96},size={149,15},bodyWidth=60,title="Incoherent XS (cm)" 875 // SetVariable MC_setvar0_2,limits={inf,inf,0.1},value= root:Packages:NIST:SAS:gSig_incoh 876 // SetVariable MC_setvar0_3,pos={491,142},size={150,15},bodyWidth=60,title="Sample Radius (cm)" 877 // SetVariable MC_setvar0_3,limits={inf,inf,0.1},value= root:Packages:NIST:SAS:gR2 878 // PopupMenu MC_popup0,pos={476,13},size={165,20},proc=MC_ModelPopMenuProc,title="Model Function" 879 // PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()" 880 // Button MC_button0,pos={480,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation" 881 // Button MC_button1,pos={644,181},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D" 882 // SetVariable setvar0_3,pos={568,484},size={50,20},disable=1 883 // GroupBox group0,pos={478,42},size={267,130},title="Monte Carlo" 884 // SetVariable cntVar,pos={653,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)" 885 // SetVariable cntVar,format="%d" 886 // SetVariable cntVar,limits={1,10,1},value= root:Packages:NIST:SAS:gCntTime 887 // 888 // String fldrSav0= GetDataFolder(1) 889 // SetDataFolder root:Packages:NIST:SAS: 890 // Edit/W=(476,217,746,450)/HOST=# results_desc,results 891 // ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150 892 // SetDataFolder fldrSav0 893 // RenameWindow #,T_results 894 // SetActiveSubwindow ## 895 EndMacro 896 897 // as a standalone panel, extra control bar (right) and subwindow implementations don't work right 898 // for various reasons... 804 899 Window MC_SASCALC() : Panel 805 900 PauseUpdate; Silent 1 // building window... 806 NewPanel /W=(92,556, 390,1028)/K=1 as "SANS Simulator"901 NewPanel /W=(92,556,713,818)/K=1 as "SANS Simulator" 807 902 SetVariable MC_setvar0,pos={28,73},size={144,15},bodyWidth=80,title="# of neutrons" 808 903 SetVariable MC_setvar0,format="%5.4g" … … 817 912 PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()" 818 913 Button MC_button0,pos={17,181},size={130,20},proc=MC_DoItButtonProc,title="Do MC Simulation" 819 Button MC_button1,pos={181,181},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D" 914 Button MC_button0,fColor=(3,52428,1) 915 Button MC_button1,pos={17,208},size={80,20},proc=MC_Display2DButtonProc,title="Show 2D" 820 916 SetVariable setvar0_3,pos={105,484},size={50,20},disable=1 821 917 GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo" 822 918 SetVariable cntVar,pos={190,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)" 823 919 SetVariable cntVar,format="%d" 824 SetVariable cntVar,limits={1,10,1},value= root:Packages:NIST:SAS:gCntTime 920 SetVariable cntVar,limits={1,60,1},value= root:Packages:NIST:SAS:gCntTime 921 Button MC_button2,pos={17,234},size={100,20},proc=SaveAsVAXButtonProc,title="Save 2D VAX" 922 CheckBox check0,pos={216,180},size={68,14},title="Raw counts",variable = root:Packages:NIST:SAS:gRawCounts 923 CheckBox check0_1,pos={216,199},size={60,14},title="Yes Offset",variable= root:Packages:NIST:SAS:gDoTraceOffset 825 924 826 925 String fldrSav0= GetDataFolder(1) 827 926 SetDataFolder root:Packages:NIST:SAS: 828 Edit/W=( 13,217,283,450)/HOST=# results_desc,results927 Edit/W=(344,23,606,248)/HOST=# results_desc,results 829 928 ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150 830 929 SetDataFolder fldrSav0 … … 832 931 SetActiveSubwindow ## 833 932 EndMacro 834 835 933 836 934 Function CountTimeSetVarProc(sva) : SetVariableControl … … 956 1054 End 957 1055 1056 // writes out a VAX binary data file 1057 // automatically generates a name 1058 // will prompt for the sample label 1059 // 1060 Function SaveAsVAXButtonProc(ctrlName) : ButtonControl 1061 String ctrlName 1062 1063 WriteVAXData("SAS","",0) 1064 End 1065 1066 958 1067 959 1068
Note: See TracChangeset
for help on using the changeset viewer.