source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/Polarization/Pol_PolCorr.ipf @ 838

Last change on this file since 838 was 816, checked in by srkline, 12 years ago

in 2D models, adjusted the point ranges to be integers. non-integer values occasionally causes issues.

changed names of polarization procedures to start with "Pol_" so that they are easier to locate.

added MixedDumbbell_FFT model as another example of a FFT fitting function.
File size: 15.9 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3// this is version 3.0? of the code for the DLL from K. Krycka
4// 18 MAY 2011
5//
6// ocnverted to Igor, SRK
7//
8//
9
10
11// this is one monolithic program to gather the inputs and then do the corrections
12//
13// -- some changes here --
14//
15// (1) data here is assumed to be RAW data, not partially corrected ASCII
16//
17// look for HARD WIRED VALUES, and eventually get rid of these
18//
19
20
21Function PolCorr_main()
22
23        Variable Monitor_Normalization = 1.0e+08 
24        Variable Max_Data_Multiplicity = 5
25        Variable Cells = 3
26        Variable Number_Cells = Cells
27        Variable Orientation = 4
28        Variable Q_PixelsX = 128
29        Variable Q_PixelsY = 128
30       
31//      Variable lenFile = 100
32//      String File_PP1,File_PP2,File_PP3,File_PP4,File_PP5
33//      String File_MP1,File_MP2,File_MP3,File_MP4,File_MP5
34//      String File_MM1,File_MM2,File_MM3,File_MM4,File_MM5
35//      String File_PM1,File_PM2,File_PM3,File_PM4,File_PM5
36//      String File_PPout,File_MPout,File_MMout,File_PMout
37        // these differed only in the capitialzation of out -> Out, which is not recognized in Igor
38        // but were apparently never used
39        //String File_PP_Out,File_MP_Out,File_MM_Out,File_PM_Out
40       
41        Variable input
42        String holder
43        Variable counter, counter2
44        Variable start_day, start_month, start_year
45        Variable helium_day, helium_month, helium_year
46        Variable data_day, data_month, data_year
47        Make/O/D/N=(Cells) He_Pol_Measure_Time, gamm, Te, onl, Initial_He_Pol
48        Make/O/D/N=(Orientation) Number_Files            //PP, MP, MM, PM
49        Make/O/T/N=(Orientation, Max_Data_Multiplicity) Data_Files
50        Make/O/T/N=(Orientation) Out_Files
51        Make/O/D/N=(Orientation,Max_Data_Multiplicity) Helium_Cell
52        Variable Polarization_SuperMirror, P_M           //polarization of supermirror
53        Variable Polarization_Flipper, P_F              //polarization of flipper = 2*efficiency - 1
54       
55        Variable YesNo
56        String ch
57        Make/O/T/N=(4) sub_holder
58        Variable ori
59        Variable mult
60        Variable Row, Column
61        String line1,line2,line3,line4,line5,line6,line7,line8,line9,line10
62        String line11,line12,line13,line14,line15,line16,line17,line18,line19
63        Variable date_v, month, year, hour, minute, second
64        Variable time_v
65        Variable moncounts
66        Variable sum_of_counts
67        Variable Scaled_Counts
68//      char *JANUARY={"JAN"};char *FEBRUARY={"FEB"};char *MARCH={"MAR"};char *APRIL={"APR"}
69//      char *MAY={"MAY"};char *JUNE={"JUN"};char *JULY={"JUL"};char *AUGUST={"AUG"}
70//      char *SEPTEMBER={"SEP"};char *OCTOBER={"OCT"};char *NOVEMBER={"NOV"};char *DECEMBER={"DEC"}
71        Variable Multiplicity
72        Variable a, b
73        Variable X, Y
74       
75        Make/O/D/N=(Orientation, Max_Data_Multiplicity) Time_Data_Collected
76        Make/O/D/N=(Orientation, Max_Data_Multiplicity) Monitor_Counts                  //[PP,MP,MM,PM][Data_Multiplicity]
77        Make/O/D/N=(Orientation, Max_Data_Multiplicity) Total_Counts                    //[PP,MP,MM,PM][Data_Multiplicity]
78        Make/O/D/N=(Orientation, Max_Data_Multiplicity, 2) He_Transmission              //[PP,MP,MM,PM][Data_Multiplicity][majority, minority] in percent transmitted.
79        Make/O/D/N=(Orientation, Q_PixelsX, Q_PixelsY) Scattering_Intensity             //[PP,MP,MM,PM][Max_Data_Multiplicity][X][Y]
80        Make/O/D/N=(Orientation, Orientation) Prefactors
81        Make/O/D/N=(Orientation, Orientation) Inverted_Matrix
82        Make/O/D/N=(Orientation) Norm_Data
83
84//***********************************************************************************************************
85//Read in Record.txt parameters
86//***********************************************************************************************************
87
88   
89////// this looks like all of the necessary inputs to make it work
90    //
91//    fscanf(user_input, "%f", &input); start_day = input; //printf("start day is %u\n", start_day)
92//    fscanf(user_input, "%f", &input); start_month = input; //printf("start month is %u\n", start_month)
93//    fscanf(user_input, "%f", &input); start_year = input; //printf("start year is %u\n", start_year)
94   
95//// HARD WIRED VALUES
96        start_day = 26
97        start_month = 4
98        start_year = 2010
99   
100    helium_day = start_day
101    helium_year = start_year - 2000
102    if(start_month == 1)
103        helium_month = 0
104    endif
105    if(start_month == 2)
106        helium_month = 31
107    endif
108    if(start_month == 3)
109            helium_month = 59
110    endif
111    if(start_month == 4)
112            helium_month = 90
113         endif
114    if(start_month == 5)
115            helium_month = 120
116         endif
117    if(start_month == 6)
118            helium_month = 151
119         endif
120    if(start_month == 7)
121            helium_month = 181
122         endif
123    if(start_month == 8)
124            helium_month = 212
125         endif
126    if(start_month == 9)
127            helium_month = 243
128         endif
129    if(start_month == 10)
130            helium_month = 273
131         endif
132    if(start_month == 11)
133            helium_month = 304
134         endif
135    if(start_month == 12)
136            helium_month = 334
137         endif
138
139         
140//// HARD WIRED VALUES
141        // right now, Cells = 1, only one cell in use
142        //
143        Cells = 1
144       
145       
146    for(counter=0; counter<Cells; counter+=1)
147        He_Pol_Measure_Time[counter] = 14.5                     //hours, HARD WIRED (from midnight of start day?)
148    endfor
149   
150    for(counter=0; counter<Cells; counter+=1)
151        gamm[counter] = 298.5           //hours
152    endfor
153   
154    for(counter=0; counter<Cells; counter+=1)
155        Te[counter] = 0.87
156    endfor
157   
158    // onl == mu == opacity of the cell
159    for(counter=0; counter<Cells; counter+=1)
160        onl[counter] = 3.108     
161    endfor
162   
163    for(counter=0; counter<Cells; counter+=1)
164        Initial_He_Pol[counter] = 0.641
165    endfor
166   
167// number of input files in each orientation
168    Number_Files[0] = 1
169    Number_Files[1] = 2
170    Number_Files[2] = 1
171    Number_Files[3] = 2
172   
173// index to keep track of which cell was used --for which data set-- at which orientation
174        Helium_Cell = 0
175//      Helium_Cell[orientation][data set num] = cell number
176        Helium_Cell[0][0] = 1
177        Helium_Cell[1][0] = 1
178        Helium_Cell[1][1] = 1
179        Helium_Cell[2][0] = 1
180        Helium_Cell[3][0] = 1
181        Helium_Cell[3][1] = 1
182 
183//     for(counter=0; counter<Orientation; counter++){
184//      for(counter2=0; counter2<Max_Data_Multiplicity; counter2++){
185//              fscanf(user_input, "%f", &input); Helium_Cell[counter][counter2] = input; //printf("  cell used %u %u is %f", counter, counter2, Helium_Cell[counter][counter2]);             
186//      }
187//    }
188   
189//    fscanf(user_input, "%f", &input); Polarization_SuperMirror = input; P_M = input; //printf("polarization of super mirror is %f \n", Polarization_SuperMirror)
190//    fscanf(user_input, "%f", &input); Polarization_Flipper = 2.0*input -1.0; P_F = input;  //printf("polarization of flipper is %f \n", Polarization_Flipper)
191
192//// HARD WIRED VALUES
193        P_M = 0.937
194        Polarization_SuperMirror = P_M
195
196        P_F = 0.98
197        Polarization_Flipper = 2*P_F -1.0
198
199
200//***********************************************************************************************************
201//Read in data files, add scans, calculate He transmissions
202//***********************************************************************************************************
203
204//
205// -- use AddFilesInList() here to put the data in SAM,
206// then move each orienatation into the Scattering_Intensity[][][] matrix
207// and fill in the He cell timing information, based on the header information.
208//
209// unroll the ori loop, and do each separtely, for better control of which files to load
210// - presumably input from a panel somewhere as run nubmers.
211//
212//      samStr = ParseRunNumberList("1,2,3,")           //returns a string with the file names
213//  err = AddFilesInList(activeType,samStr)                     // adds the files in samStr to the activeType folder
214//
215// String fname = FindFileFromRunNumber(num) // could also be used. this returns the full path:name
216//
217//
218// -- I'm still not exactly sure what the times really are and how they are calculated. write this up...
219//
220//
221//      PathInfo catPathName                    //this is where the files are
222//      String pathStr=S_path
223               
224    Scattering_Intensity = 0            //end row and column initialization loops
225
226        for(ori = 0; ori<4; ori+=1)                             //[PP,MP,MM,PM]
227                Multiplicity=Number_Files[ori]          //[SS,ST,CS,CT][PP,MP,MM,PM]
228               
229                for(mult = 0; mult<Multiplicity; mult+=1)               //Multiplicity = number of relevant files   
230       
231                       
232                        Monitor_Counts[ori][mult] = moncounts
233                       
234
235// for this time value, get the midpoint of the data set, that is start + 1/2 of the collection time
236                        Time_Data_Collected[ori][mult] = time_v
237
238                       
239                        Variable He_cell = Helium_Cell[ori][mult]
240                       
241                        Variable He_pol_time = Initial_He_Pol[He_cell]*exp(-(Time_Data_Collected[ori][mult]-He_Pol_Measure_Time[He_cell])/gamm[He_cell])
242                        He_Transmission[ori][mult][0] = Te[He_cell]*exp(-(1-He_pol_time)*onl[He_cell])       //He majority transmission
243                        He_Transmission[ori][mult][1] = Te[He_cell]*exp(-(1+He_pol_time)*onl[He_cell])   //He minority transmission
244                       
245                        sum_of_counts = 0 
246                        for(Row=0; Row<128; Row+=1)
247                                for(Column=0; Column<128; Column+=1)                   
248                                        Scaled_Counts = input*Monitor_Normalization/Monitor_Counts[ori][mult]           //[PP,MP,MM,PM][Data_Multiplicity][X][Y]
249                                        Scattering_Intensity[ori][Column][Row] = Scattering_Intensity[ori][Column][Row] + Scaled_Counts
250                                        sum_of_counts = sum_of_counts + Scaled_Counts 
251                                endfor
252                        endfor
253                        Total_Counts[ori][mult] = sum_of_counts
254                               
255                endfor          //end Multiplicity loop
256        endfor                  //end Orientation loop
257       
258       
259
260//***********************************************************************************************************
261//Fill in matrix prefactors prior to solving four simultaneous equations
262//***********************************************************************************************************
263
264
265//
266// should this be a time-of-collection weighted linear combination of the coefficients?
267// what if each data file that contributes to Iab(q) was not collected for the same length
268// of time?
269//
270//
271//
272//
273        for(a=0; a<4; a+=1)
274                for(b=0; b<4; b+=1)
275                        Prefactors[a][b] = 0.0
276                endfor
277        endfor
278
279        Multiplicity = Number_Files[0]                   //[PP]
280        for(mult = 0; mult<Multiplicity; mult+=1)
281                Prefactors[0][0] = Prefactors[0][0] + (1+P_M)*He_Transmission[0][mult][0]
282                Prefactors[0][1] = Prefactors[0][1] + (1-P_M)*He_Transmission[0][mult][0]
283                Prefactors[0][2] = Prefactors[0][2] + (1-P_M)*He_Transmission[0][mult][1]
284                Prefactors[0][3] = Prefactors[0][3] + (1+P_M)*He_Transmission[0][mult][1]
285        endfor
286
287        Multiplicity = Number_Files[1]                   //[MP]
288        for(mult = 0; mult<Multiplicity; mult+=1)
289                Prefactors[1][0] = Prefactors[1][0] + (1-P_M*P_F)*He_Transmission[1][mult][0]
290                Prefactors[1][1] = Prefactors[1][1] + (1+P_M*P_F)*He_Transmission[1][mult][0]
291                Prefactors[1][2] = Prefactors[1][2] + (1+P_M*P_F)*He_Transmission[1][mult][1]
292                Prefactors[1][3] = Prefactors[1][3] + (1-P_M*P_F)*He_Transmission[1][mult][1]
293        endfor
294
295        Multiplicity = Number_Files[2]                  //[MM]
296        for(mult = 0; mult<Multiplicity; mult+=1)
297                Prefactors[2][0] = Prefactors[2][0] + (1-P_M*P_F)*He_Transmission[2][mult][1]
298                Prefactors[2][1] = Prefactors[2][1] + (1+P_M*P_F)*He_Transmission[2][mult][1]
299                Prefactors[2][2] = Prefactors[2][2] + (1+P_M*P_F)*He_Transmission[2][mult][0]
300                Prefactors[2][3] = Prefactors[2][3] + (1-P_M*P_F)*He_Transmission[2][mult][0]
301        endfor
302
303        Multiplicity = Number_Files[3]                  //[PM]
304        for(mult = 0; mult<Multiplicity; mult+=1)
305                Prefactors[3][0] = Prefactors[3][0] + (1+P_M)*He_Transmission[3][mult][1]
306                Prefactors[3][1] = Prefactors[3][1] + (1-P_M)*He_Transmission[3][mult][1]
307                Prefactors[3][2] = Prefactors[3][2] + (1-P_M)*He_Transmission[3][mult][0]
308                Prefactors[3][3] = Prefactors[3][3] + (1+P_M)*He_Transmission[3][mult][0]
309        endfor
310
311//***********************************************************************************************************
312//Written by K. Krycka 2006, modified for SANS Nov. 2007.
313//This is an n-dimensional square matrix solver.
314//***********************************************************************************************************
315
316        Variable subtract
317        Variable col, inner, inner_col, inner_row
318       
319        Variable dimension = 4
320        Make/O/D/N=(4,4) Identity
321        Identity = 0
322        Identity[0][0] = 1
323        Identity[1][1] = 1
324        Identity[2][2] = 1
325        Identity[3][3] = 1
326       
327        Make/O/D/N=(dimension, dimension) MatrixToInvert
328        for(col = 0; col<dimension; col+=1)
329                for(row = 0; row<dimension; row+=1)
330                MatrixToInvert[row][col] = Prefactors[row][col]
331                        if(col==row)
332                                Inverted_Matrix[row][col] = 1.0
333                        else
334                                Inverted_Matrix[row][col] = 0.0
335                        endif
336                endfor
337        endfor
338       
339        subtract = 0
340        if(MatrixToInvert[0][0]!=0)
341                for(col = 0; col<dimension-1; col+=1)
342                        for(row = col+1; row<dimension; row+=1)
343                                subtract = MatrixToInvert[row][col] / MatrixToInvert[col][col]
344                                if(MatrixToInvert[row][col]!=0)
345                                        for(inner_col=0; inner_col<dimension; inner_col+=1)
346                                                Identity[row][inner_col] = Identity[row][inner_col] - subtract*Identity[col][inner_col]
347                                                MatrixToInvert[row][inner_col] = MatrixToInvert[row][inner_col] - subtract*MatrixToInvert[col][inner_col]
348                                        endfor
349                                endif
350                        endfor
351                endfor
352       
353                       
354                for(row = 0; row<dimension-1; row+=1)
355                        for(col = row+1; col<dimension; col+=1)
356                                subtract = MatrixToInvert[row][col] / MatrixToInvert[col][col]
357                                if(MatrixToInvert[row][col]!=0)
358                                        for(inner_row=0; inner_row<dimension; inner_row+=1)
359                                                Identity[row][inner_row] = Identity[row][inner_row] - subtract*Identity[col][inner_row]
360                                                MatrixToInvert[row][inner_row] = MatrixToInvert[row][inner_row] - subtract*MatrixToInvert[col][inner_row]
361                                        endfor
362                                endif
363                        endfor
364                endfor
365       
366        endif
367       
368        Make/O/D/N=(dimension) divide
369    for(row=0; row<dimension; row+=1)
370        divide[row] =  MatrixToInvert[row][row]
371    endfor
372
373    for(row=0; row<dimension; row+=1)
374        for(col=0; col<dimension; col+=1)
375                Identity[row][col] = Identity[row][col]/divide[row]
376                MatrixToInvert[row][col] = MatrixToInvert[row][col]/divide[row]
377        endfor
378    endfor
379
380        for(a=0; a<dimension; a+=1)
381                for(b=0; b<dimension; b+=1)
382                        Inverted_Matrix[a][b] = Identity[a][b]
383                endfor
384        endfor
385
386//***********************************************************************************************************
387//Solve for cross sections PP, MP, MM, and PM
388//***********************************************************************************************************
389//
390// it looks like the Scattering_Intensity matrix is overwritten in place
391//
392        for(X=0; X<Q_PixelsX; X+=1)
393                for(Y=0; Y<Q_PixelsY; Y+=1)
394               
395                        for(a=0; a<4; a+=1)
396                                Norm_Data[a] = Scattering_Intensity[a][X][Y]
397                        endfor
398                       
399                        for(a=0; a<4; a+=1)
400                                Scattering_Intensity[a][X][Y]=0
401                        endfor
402               
403                        for(ori=0; ori<Orientation; ori+=1)
404                                for(col=0; col<4; col+=1)
405                               
406                                        Scattering_Intensity[ori][X][Y] = Scattering_Intensity[ori][X][Y] + Inverted_Matrix[ori][col]*Norm_Data[col]
407                                        //Solved[ori][X][Y] = Solved[ori][X][Y] + Inverted_Matrix[ori][col]*Scattering_Intensity[col][X][Y]
408                               
409                                endfor                  //end col loop
410                        endfor                  //end ori loop
411               
412                endfor                  //end Y loop
413        endfor                  //end X loop
414
415//***********************************************************************************************************
416//Write processed, output files
417//***********************************************************************************************************
418
419
420// CHANGE this step to write out the 4 cross sections from the Scattering_Intensity[][][] matrix, just to simple
421// arrays, to be carried through the reduction
422
423//
424        Make/O/D/N=(Q_PixelsX,Q_PixelsY) linear_data_PP, linear_data_MP, linear_data_MM, linear_data_PM
425       
426        linear_data_PP = Scattering_Intensity[0][p][q]
427        linear_data_MP = Scattering_Intensity[1][p][q]
428        linear_data_MM = Scattering_Intensity[2][p][q]
429        linear_data_PM = Scattering_Intensity[3][p][q]
430
431       
432//      for(ori=0; ori<Orientation; ori+=1)
433//     
434//              for(Row = 0; Row<Q_PixelsX; Row+=1)
435//                      for(Column = 0; Column<Q_PixelsY; Column+=1)
436//                              //fprintf(fp, "%f\n", Raw_Counts[ori][0][Row][Column])
437//                              //fprintf(fp, "%f\n", Scattering_Intensity[ori][Row][Column])
438//                              fprintf(fp, "%f\n", Scattering_Intensity[ori][Row][Column])
439//                      endfor
440//              endfor
441//             
442//      endfor          //end ori loop
443//
444//
445//
446
447
448//***********************************************************************************************************
449//Write message to user of completion
450//***********************************************************************************************************
451   
452
453    print "start year", start_year
454    print "Time of first data collection is ", Time_Data_Collected[0][0]
455    print "Monitor counts are ", Monitor_Counts[0][0]
456    print "Total counts are ", Total_Counts[0][0]
457   
458        Print "Your data has beed processed"   
459
460End  // PolCorr_main
Note: See TracBrowser for help on using the repository browser.