Ignore:
Timestamp:
Feb 3, 2009 3:57:08 PM (14 years ago)
Author:
srkline
Message:

A large number of changes and fixes:

--168/169/170: panels and windows are now at least on-screen for all packages. Tested
with 1024x768 resolution.
-- closed ticket 176 which was a question about resampling data to generate error estimates
on fitted parameters. Useful for reflectometry, not needed for SANS.
--157: bug of low Q power law extrapolation in Invariant fixed by avoiding q==0
--178/180: Tr/Tw? notification in USANS. log/lin checkbox for display.
--167: saveData checkbox for USANS not behaving well. turns off/on better now.
--197: changed all (?) 1D writing routines to enforce 26 characters as the maximum length
to make sure that file loading will never exceed 31 characters

-- lots of changes to MonteCarlo? and SASCALC

  • SASCALC now enforces *exact* lens conditions, rather than a close approximation
  • improved MonteCarlo? interface panel
  • added writer for simlulated VAX binary data file
  • can save 2D as ABS or raw counts
  • can freeze w/no offset
File:
1 edited

Legend:

Unmodified
Added
Removed

  • sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/MultScatter_MonteCarlo_2D.ipf

    r465 r472  
    5656// - a number of models are now found to be ill-behaved when q=1e-10. Then the random deviate calculation blows up. 
    5757//   a warning has been added - but some models need a proper limiting value, and some (power-law) are simply unuseable 
    58 //   unless something else can be done. 
     58//   unless something else can be done. Using a log-spacing of points doesn't seem to help, and it introduces a lot of 
     59//   other problems. Not the way to go. 
    5960// - if the MC gags on a simulation, it often gets "stuck" and can't do the normal calculation from the model, which it 
    6061//   should always default to... 
    … …  
    7071        // does nothing in the Igor code 
    7172        Duplicate/O results retWave 
    72         //results[0] = -1*(datetime) 
    7373 
    7474        Variable NNeutron=inputWave[0] 
    … …  
    8181         
    8282        variable mt= ThreadGroupCreate(nthreads) 
     83        NVAR gInitTime = root:Packages:NIST:SAS:gRanDateTime            //time that SASCALC was started 
     84 
    8385         
    8486        inputWave[0] = NNeutron/nthreads                //split up the number of neutrons 
    … …  
    9395                Duplicate/O inputWave $("inputWave"+num2istr(i)) 
    9496                Duplicate/O ran_dev $("ran_dev"+num2istr(i)) 
    95                  
     97                                 
    9698                // ?? I need explicit wave references? 
    9799                // maybe I need to have everything in separate data folders - bu tI haven't tried that. seems like a reach. 
    … …  
    99101                if(i==0) 
    100102                        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*(datetime-gInitTime)           //to initialize ran3 
    102105                        ThreadStart mt,i,Monte_SANS_W1(inputWave0,ran_dev0,nt0,j10,j20,nn0,linear_data0,retWave0) 
    103106                        Print "started thread 0" 
    … …  
    105108                if(i==1) 
    106109                        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*(datetime-gInitTime)           //to initialize ran1 
    108112                        ThreadStart mt,i,Monte_SANS_W2(inputWave1,ran_dev1,nt1,j11,j21,nn1,linear_data1,retWave1) 
    109113                        Print "started thread 1" 
    … …  
    180184 
    181185// worker function for threads, does nothing except switch between XOP and Igor versions 
     186// 
     187// uses ran3 
    182188ThreadSafe Function Monte_SANS_W1(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
    183189        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 
    … …  
    191197        return (0) 
    192198End 
     199 
    193200// worker function for threads, does nothing except switch between XOP and Igor versions 
     201// 
     202// uses ran1 
    194203ThreadSafe Function Monte_SANS_W2(inputWave,ran_dev,nt,j1,j2,nn,linear_data,results) 
    195204        WAVE inputWave,ran_dev,nt,j1,j2,nn,linear_data,results 
    … …  
    569578        qu = 4*pi/lam            
    570579         
    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:gQ 
    573 //      WAVE xw = root:Packages:NIST:SAS:xw 
    574  
    575580// hard-wired into the Simulation directory rather than the SAS folder. 
    576581// plotting resolution-smeared models won't work any other way 
    … …  
    578583        WAVE Gq = root:Simulation:gQ 
    579584        WAVE xw = root:Simulation:xw 
    580         SetScale/I x (0+1e-6),qu*(1-1e-10),"", 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+1e-4),qu*(1-1e-10),"", Gq,xw                    //don't start at zero or run up all the way to qu to avoid numerical errors 
     586 
    583587/// if all of the coefficients are well-behaved, then the last point is the background 
    584588// and I can set it to zero here (only for the calculation) 
    … …  
    593597        Gq = Gq*sin(2*asin(x/qu))/sqrt(1-(x/qu))                        // exact 
    594598        // 
     599        // 
    595600        Integrate/METH=1 Gq/D=Gq_INT 
    596601         
    … …  
    605610End 
    606611 
    607  
     612// returns the random deviate as a wave 
     613// and the total SAS cross-section [1/cm] sig_sas 
     614// 
     615// uses a log spacing of x for better coverage 
     616// downside is that it doesn't use built-in 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 power-laws) 
     622// giving unreasonably large SAS cross sections. (>>10) 
     623// 
     624Function 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=1e-5 
     633        qu = 4*pi/lam            
     634 
     635// hard-wired into the Simulation directory rather than the SAS folder. 
     636// plotting resolution-smeared 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+1e-4),qu*(1-1e-10),"", 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 well-behaved, 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[num-1] = 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_ran-1] 
     660         
     661        Gq_INT /= Gq_INT[nPts_ran-1] 
     662         
     663        Duplicate/O Gq_INT $outWave 
     664 
     665        return(0) 
     666End 
    608667 
    609668ThreadSafe Function FindPixel(testQ,testPhi,lam,sdd,pixSize,xCtr,yCtr,xPixel,yPixel) 
    … …  
    802861 
    803862// 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 ## 
     895EndMacro 
     896 
     897// as a stand-alone panel, extra control bar  (right) and subwindow implementations don't work right  
     898// for various reasons... 
    804899Window MC_SASCALC() : Panel 
    805900        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" 
    807902        SetVariable MC_setvar0,pos={28,73},size={144,15},bodyWidth=80,title="# of neutrons" 
    808903        SetVariable MC_setvar0,format="%5.4g" 
    … …  
    817912        PopupMenu MC_popup0,mode=1,value= #"MC_FunctionPopupList()" 
    818913        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" 
    820916        SetVariable setvar0_3,pos={105,484},size={50,20},disable=1 
    821917        GroupBox group0,pos={15,42},size={267,130},title="Monte Carlo" 
    822918        SetVariable cntVar,pos={190,73},size={80,15},proc=CountTimeSetVarProc,title="time(s)" 
    823919        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 
    825924         
    826925        String fldrSav0= GetDataFolder(1) 
    827926        SetDataFolder root:Packages:NIST:SAS: 
    828         Edit/W=(13,217,283,450)/HOST=#  results_desc,results 
     927        Edit/W=(344,23,606,248)/HOST=#  results_desc,results 
    829928        ModifyTable format(Point)=1,width(Point)=0,width(results_desc)=150 
    830929        SetDataFolder fldrSav0 
    … …  
    832931        SetActiveSubwindow ## 
    833932EndMacro 
    834  
    835933 
    836934Function CountTimeSetVarProc(sva) : SetVariableControl 
    … …  
    9561054End 
    9571055 
     1056// writes out a VAX binary data file 
     1057// automatically generates a name 
     1058// will prompt for the sample label 
     1059// 
     1060Function SaveAsVAXButtonProc(ctrlName) : ButtonControl 
     1061        String ctrlName 
     1062 
     1063        WriteVAXData("SAS","",0) 
     1064End 
     1065 
     1066 
    9581067 
    9591068 
Note: See TracChangeset for help on using the changeset viewer.