pro calc_eff_area_nustar_v11, photon_file, earr, fmodel, mode, bsel=bsel, sigma=sigma, $ rsigma=rsigma, outfile=outfile, resfile=resfile, atten=atten, satten=satten, $ outdir=outdir, mon=mon, filemon=filemon, rrxorough=rrxorough, rdmin=rdmin, rdmax=rdmax ; ; Using a photon file from MT_RAYOR (mt_save,mode="p") to deduce ; NuSTAR effective area ; ; By the formula: Eff_area = sum_over_good_photons(r_1_i * r_2_i) * aperture_area / number_of_injected_photons ; common master_table, MTF_glass,MTF_location,MTF_optic, $ MTF_layer, MTF_uplo, MTF_pos_in_optic, MTF_run_number, MTF_recipe, MTF_witness, $ MTF_poschamb, MTF_separ, MTF_mount_plate, MTF_pos_on_plate if n_params() lt 4 then begin print,'Syntax: calc_eff_area_nustar_v11, photon_file, earr, fmodel, mode, sigma=, outfile=,' print,' resfile=, outdir=, atten=, mon=, rsigma=, rrxorough=' print,'' print,' earr [keV] is energy array.' print,' fmodel is flight model number (1 or 2).' print,' mode tells what to include in reflectivity calculations (see below).' print,' Keywords' print,' ''bsel'' will cause a selection on the bounce value print,' ''sigma'' [A] will replace the roughness value except for' print,' ''as-coated'' when witness fits are available.' print,' ''rsigma'' a factor to apply to all roughnesses.' print,' ''rrxorough'' a factor to apply to all RXO coated roughnesses.' print,' ''outfile'' will add this to the output filename.' print,' ''resfile'' will make this the output filename (overrides outfile).' print,' ''outdir'' will cause outfiles to be written there.' print,' ''atten'' will only use one out of ''atten'' photons.' print,' ''satten'' is the subselection under ''atten''' print,' 0 <= satten < atten' print,' ''mon'' is for producing a monitor file print,' ''rdmin'' relative change factor for all ''as-coated'' recipes' print,' ''rdmax'' relative change factor for all ''as-coated'' recipes' print,'' print,' mode must be a four character string with 0s or 1s as' print,' ''rgda'' signifying :' print,' ''a'' as-coated (1) or standard recipes (0)' print,' ''d'' include d-spacing variations (1) or disregard (0)' print,' ''g'' include gamma variations (1) or disregard (0)' print,' ''r'' include roughness variations (1) or disregard (0)' print,' Example: ''0101'' will use witness sample fits and just gamma variations' print,'' print,' Only status zero photons will be included and an additional selection' print,' on the bounce value may be imposed by the ''bsel'' keyword.' print,'' return endif ; -- Version number (string) version = '11' if fmodel lt 1 or fmodel gt 2 then begin print,'The ''fmodel'' argument was outside allowed range' return endif mstr = 'FM'+stri(fmodel) sz = size(mode) if sz(0) ne 0 or sz(1) ne 7 then begin print,'''mode'' is not a scalar string' return endif if not is_digit(mode) then begin print,'''mode'' holds non-digits' return endif use_as_coated = long(strmid(mode,3,1)) use_dmod_var = long(strmid(mode,2,1)) use_gamma_var = long(strmid(mode,1,1)) use_rough_var = long(strmid(mode,0,1)) if ~keyword_set(rsigma) then rsigma = 1.0 if ~keyword_set(rrxorough) then rrxorough = 1.0 if ~keyword_set(rdmin) then rdmin = 1.0 if ~keyword_set(rdmax) then rdmax = 1.0 if keyword_set(mon) then mon = 1 else mon = 0 ; set the attenuation and satten if not keyword_set(atten) then atten = 1 if keyword_set(satten) then begin if satten ge atten then satten = atten - 1 if satten lt 0 then satten = 0 endif else begin satten = 0 endelse astr = string(atten,satten,format='("_at",i03,"i",i03)') ; define output filename if keyword_set(outdir) then begin outdir = app_slash(outdir) endif else begin outdir = '' endelse if keyword_set(outfile) then begin outfilex = '_'+outfile endif else outfilex = '' if keyword_set(bsel) then begin bstr = '_b'+string(bsel,format='(i1)') endif else bstr = '' if keyword_set(resfile) then begin outname = outdir+resfile endif else begin outname = outdir+'eff_area_calc_v'+version+'_'+mstr+'_mode'+mode+astr+outfilex+bstr+'.fits' endelse if mon then begin if keyword_set(filemon) then begin filemon = outdir+filemon endif else begin filemon = outdir+'monitor_v'+version+'_'+mstr+'_mode'+mode+astr+outfilex+bstr+'.fits' endelse endif ; get the basic information from the photon file dol = photon_file+'+1' pfhdr = headfits( photon_file, ext=1 ) aperture_area = fxpar( pfhdr, 'fraparea' ) n_inject = long(fxpar( pfhdr, 'n_inject' )) sampling = fxpar( pfhdr, 'sampling' ); Take the original sampling or attenuation n_inject = n_inject / sampling ; into account ; read the photon properties from the photon file rdfitscol, dol, mirror_arr, 'mirror', nrows=nrows, /silent rdfitscol, dol, ang1_arr, 'angle_in1', /silent rdfitscol, dol, ang2_arr, 'angle_in2', /silent rdfitscol, dol, i1z_arr, 'i1z', /silent; z value of impact position on mirror rdfitscol, dol, i2z_arr, 'i2z', /silent rdfitscol, dol, azimuth_arr, 'azimuth', /silent rdfitscol, dol, status_arr, 'status', /silent rdfitscol, dol, bounce_arr, 'bounce', /silent w = where( status_arr eq 0, nw ) if nw eq 0 then begin print,'No status zero photons found ...' return endif mirror_arr = mirror_arr(w) ang1_arr = ang1_arr(w) ang2_arr = ang2_arr(w) i1z_arr = i1z_arr(w) i2z_arr = i2z_arr(w) azimuth_arr = azimuth_arr(w) bounce_arr = bounce_arr(w) if keyword_set(bsel) then begin w = where( bounce_arr eq bsel, nw ) if nw eq 0 then begin print,'No photons found with specified bounce value ...' return endif mirror_arr = mirror_arr(w) ang1_arr = ang1_arr(w) ang2_arr = ang2_arr(w) i1z_arr = i1z_arr(w) i2z_arr = i2z_arr(w) azimuth_arr = azimuth_arr(w) bounce_arr = bounce_arr(w) endif nrows = nw; The reduced number of photons ; subtract 150 deg from azimuth to conform to the NuSTAR definition azimuth_arr = zero2pi(azimuth_arr - 150.*!pi/180.) aperture_area = aperture_area * 0.01 ; convert from mm2 to cm2 ; Define arrays to receive the monitor output if mon then begin n_out_phots = (nrows + atten - 1 - satten)/atten mon_glass = intarr( n_out_phots ) mon_layer = intarr( n_out_phots ) mon_witness_u = strarr( n_out_phots ) mon_witness_l = strarr( n_out_phots ) ; -- special for extra debugging BEGIN ; find 7 indices for the energies 10, 20, ..., 70 mon_idx = intarr(7) for ie = 1, 7 do begin w = where( abs(earr - ie*10) eq min(abs(earr - ie*10)) ) mon_idx(ie-1) = w(0) endfor mon_s_num = intarr( n_out_phots); sextant number (sector - 1) mon_t_num = intarr( n_out_phots); twelftant number (subsector - 1) mon_i1z = fltarr( n_out_phots); z of impact in 'U' mon_i2z = fltarr( n_out_phots); z of impact in 'L' mon_ang1 = fltarr( n_out_phots); grazing angle of impact in 'U' mon_ang2 = fltarr( n_out_phots); grazing angle of impact in 'L' mon_rco1 = fltarr( 7, n_out_phots); coefficient of reflection in 'U' mon_rco2 = fltarr( 7, n_out_phots); coefficient of reflection in 'L' mon_z1_rel = fltarr( n_out_phots); relative z of impact in 'U' mon_z2_rel = fltarr( n_out_phots); relative z of impact in 'L' mon_phi = fltarr( n_out_phots); Photon azimuth mon_phi_rel6 = fltarr( n_out_phots); Relative photon azimuth in sectors mon_phi_rel12 = fltarr( n_out_phots); Relative photon azimuth in subsectors mon_dmod1 = fltarr( n_out_phots ); d-spacing modification mon_dmod2 = fltarr( n_out_phots ); d-spacing modification mon_mtfu = intarr( n_out_phots ); Mirror index in Master Table, upper mirror mon_zhat1 = fltarr( n_out_phots ); zhat1 (dmods x-value) mon_zhat2 = fltarr( n_out_phots ); zhat2 (dmods x-value) mon_phihat6 = fltarr( n_out_phots ); phihat6 (dmods phi-value) mon_phihat12 = fltarr( n_out_phots ); phihat12 (dmods phi-value) ; ; All flags are zero if OK. ; flag_catdrop = 1 if there has been one or more cathode dropouts ; flag_baddata = 0 if chi2cut <= 8, else the chi2cut value ; flag_badfit = 0 if chi2 <= 1.2, else 10*chi2 value ; ; To cover both reflections the 'U' or 1 (one) reflection flag value ; is multiplied with 200 (this value because largest flag value is 139 ; and we must stay inside int range up to 32767. The 'L' or 2 (two) ; flag value is simply added. ; mon_flag_catdrop = intarr( n_out_phots ); Flag for cathode dropouts mon_flag_baddata = intarr( n_out_phots ); Flag for poor data mon_flag_badfit = intarr( n_out_phots ); Flag for poor fit mon_flag_rxo = intarr( n_out_phots ); Flag for RXO coating mon_flag_usedfit = intarr( n_out_phots ); Flag for used witness data fit ; -- special for extra debugging END endif ;--------------------------------------------------------- ; define the standard recipes recipe_std = recipestruct(); An NB function groups_start = [1,13,25,37,50,63,77,90,105,119] ;-------------- Get information from MasterTableFinal --------- read_master_table,fm=fmodel; Fills the common block: master_table ; with appropriate results ;------------------------------------------------------------------------- ; start preparations for use_as_coated == 1 ;------------------------------------------------------------------------- ;------------- read info for NB fitted witness samples ---------- ; ; Where Nicolai's fit results can be found ;;+ subdir = 'NB_data_110926/finalfits/fitsout/' ;;+ subdir = 'NB_data_120120/finalfits/fitsout/' ;;+ subdir = 'NB_data_120321/' ; Read the composite data ;;+ dol = subdir+'composite_120321.fits+1' dol = getenv('WITNESS_FIT_PARAMS_FILE')+'+1' rdfitscol, dol, nb_sample, 'SAMPLE', /silent rdfitscol, dol, angle, 'ANGLE',/silent rdfitscol, dol, nrecipe, 'NRECIPE',/silent rdfitscol, dol, c, 'C',/silent rdfitscol, dol, n, 'N',/silent rdfitscol, dol, dmin, 'DMIN',/silent rdfitscol, dol, dmax, 'DMAX',/silent rdfitscol, dol, topgam, 'TOPGAM',/silent rdfitscol, dol, gam, 'GAM',/silent rdfitscol, dol, interface, 'INTERFACE',/silent rdfitscol, dol, rosub, 'ROSUB',/silent rdfitscol, dol, roughgrad, 'ROUGHGRAD',/silent rdfitscol, dol, dtot1, 'DTOT1',/silent rdfitscol, dol, dtot2, 'DTOT2',/silent rdfitscol, dol, dtot3, 'DTOT3',/silent rdfitscol, dol, dtot4, 'DTOT4',/silent rdfitscol, dol, dtot5, 'DTOT5',/silent rdfitscol, dol, dtot6, 'DTOT6',/silent rdfitscol, dol, gtot2, 'GTOT2',/silent rdfitscol, dol, gtot3, 'GTOT3',/silent rdfitscol, dol, gtot4, 'GTOT4',/silent rdfitscol, dol, gtot5, 'GTOT5',/silent rdfitscol, dol, rough_hz, 'ROUGH_HZ',/silent rdfitscol, dol, rough_lz, 'ROUGH_LZ',/silent rdfitscol, dol, ncathodemod, 'NCATHODEMOD',/silent rdfitscol, dol, cathodemod, 'CATHODEMOD',/silent ; Update 'dmin' and 'dmax' according to keywords 'rdmin' and 'rdmax' dmin = dmin * rdmin; dmax = dmax * rdmax; ; Read the flags from 2012-03-21 rdfitscol, dol, flag_catdrop, 'FLAG_CATDROP',/silent rdfitscol, dol, flag_baddata, 'FLAG_BADDATA',/silent rdfitscol, dol, flag_badfit, 'FLAG_BADFIT',/silent ; Turn sample names to lower case to be sure that comparisons will work ; and trim the strings nb_sample = strlowcase( strtrim(nb_sample,2) ) ; Do the same thing for the witness sample names as found in the Master Table Final MTF_witness = strlowcase(strtrim( MTF_witness, 2)) lmtf_witness = strlen( MTF_witness ) ;------------------------------------------------------------------------- ; end preparations for use_as_coated == 1 ;------------------------------------------------------------------------- ;------------------------------------------------------------------------- ; Setting the roughness parameter, used throughout the stack ;------------------------------------------------------------------------- rough = recipe_std(0).roughnessHZ; This setting is applied for use_as_coated == 0 if keyword_set(sigma) then rough = sigma ;---------------------------- ; Read ML variation parabolic fit coefficients for all three parameters ; for DTU coated mirrors ;;+ afit_t_dol = '/home/njw/0idl/imd/uniformity/afit_t_120103.fits+1' afit_t_dol = getenv('UNIFORMITY_AFIT_T_FILE')+'+1' rdfitscol,afit_t_dol, afit_t_coefs, 'afit',nrows=nafit,/silent rdfitscol,afit_t_dol, afit_t_pp, 'ppos',/silent rdfitscol,afit_t_dol, afit_t_sp, 'sepplate',/silent rdfitscol,afit_t_dol, afit_t_material, 'material',/silent ;;+ afit_g_dol = '/home/njw/0idl/imd/uniformity/afit_g_120103.fits+1' afit_g_dol = getenv('UNIFORMITY_AFIT_G_FILE')+'+1' rdfitscol,afit_g_dol, afit_g_coefs, 'afit',nrows=nafit,/silent rdfitscol,afit_g_dol, afit_g_pp, 'ppos',/silent rdfitscol,afit_g_dol, afit_g_sp, 'sepplate',/silent rdfitscol,afit_g_dol, afit_g_material, 'material',/silent ;;+ afit_r_dol = '/home/njw/0idl/imd/uniformity/afit_r_120103.fits+1' afit_r_dol = getenv('UNIFORMITY_AFIT_R_FILE')+'+1' rdfitscol,afit_r_dol, afit_r_coefs, 'afit',nrows=nafit,/silent rdfitscol,afit_r_dol, afit_r_pp, 'ppos',/silent rdfitscol,afit_r_dol, afit_r_sp, 'sepplate',/silent rdfitscol,afit_r_dol, afit_r_material, 'material',/silent ; Read ML variation parabolic fit coefficients for the dspacing parameter ; for RXO coated mirrors ; There are only 2 samples and it is uncertain how they should be distributed ;;+ afit_rxo_t_dol = '/home/njw/0idl/imd/uniformity/afit_rxo_t.fits+1' afit_rxo_t_dol = getenv('UNIFORMITY_AFIT_RXO_T_FILE')+'+1' rdfitscol,afit_rxo_t_dol, afit_rxo_t_coefs, 'afit',nrows=nafit_rxo,/silent ; Define the azimuth steps for sextant and twelftant mirrors dphi6 = (2*!pi)/6.; in radians dphi12 = dphi6 / 2.; in radians ; initialize the array for summing n_earr = n_elements(earr) rcoef_sum = replicate(0.0, n_earr) ;---------------------------------------------------------------------- ; Looping over all photons -------------------------------------- ;---------------------------------------------------------------------- n_terms = 0L; A counter for the number of rcoefs, for averaging imon = -1L for iphot = long(satten), nrows-1, atten do begin ; looping over all photons n_terms = n_terms + 1L imon = imon + 1; index for monitor output if mon then mon_layer(imon) = mirror_arr(iphot) if iphot+1 mod 100 eq 0 then print,format='(i6,"/",i6)', iphot, nrows ; determine the mirror group == recipe number ; Ups - no! This is not always true. When the Master Table is inspected ; some mirrors are mixed up so it is better to get the recipe number ; from the MTF_recipe array. On the other hand the boarder beween ; Pt/C and W/Si is sharp between groups 7 and 8 (layers 89 - 90). w = where( mirror_arr(iphot) ge groups_start, nw ) group_idx = w(nw-1) group_num = group_idx + 1; counting from 1 ; In case nominal coatings are asked for the variable 'recipe_idxu' and ; 'recipe_idxl' will never be defined, so a preliminary setting will be done ; here. It will eventually be overridden if 'as-coated' is asked for. recipe_idxu = group_idx recipe_idxl = group_idx ; --- prepare detailed impact point on mirror --- sextant_num = long(azimuth_arr(iphot) / dphi6); counting starting with zero twelftant_num = long(azimuth_arr(iphot) / dphi12); counting starting with zero ; Find relative phi value phi_rel = (phi - phi1) / (phi2 - phi1) ; within the actual piece of glass ; -- For version 04 these are no longer required since the 9-tile division ; has been abandoned from this version, but they are kept for the time ; being to minimize work and they could still have some debugging value. phi_rel6 = (azimuth_arr(iphot) - sextant_num*dphi6)/ dphi6 iphi_rel6 = long(phi_rel6 * 3.) phi_rel12 = (azimuth_arr(iphot) - twelftant_num*dphi12)/ dphi12 iphi_rel12 = long(phi_rel12 * 3.) ; phi_mid6 = sextant_num*dphi6 + dphi6/2 ; note that direction of phihat is opposite Nicolai's y coordinate phihat6 = (sextant_num*dphi6 + dphi6/2 - azimuth_arr(iphot)) * (180./!pi); in degrees phihat12 = (twelftant_num*dphi12 + dphi12/2 - azimuth_arr(iphot)) * (180./!pi); in degrees ; Find the relative z position for both reflections zhat1 = -i1z_arr(iphot); opposite Nicolai's x coordinate zhat2 = -i2z_arr(iphot); opposite Nicolai's x coordinate z1_rel = (i1z_arr(iphot) + 112.5) / 225.; A mirror has a length of 225 mm ; and the z value goes from -112.5 to 112.5 iz1_rel = long(z1_rel * 3.) z2_rel = (i2z_arr(iphot) + 112.5) / 225.; A mirror has a length of 225 mm iz2_rel = long(z2_rel * 3.) if mon then begin ; -- special for extra debugging BEGIN mon_s_num(imon) = sextant_num mon_t_num(imon) = twelftant_num mon_i1z(imon) = i1z_arr(iphot) mon_i2z(imon) = i2z_arr(iphot) mon_ang1(imon) = ang1_arr(iphot) mon_ang2(imon) = ang2_arr(iphot) mon_z1_rel(imon) = z1_rel mon_z2_rel(imon) = z2_rel mon_phi(imon) = azimuth_arr(iphot) mon_phi_rel6(imon) = phi_rel6 mon_phi_rel12(imon) = phi_rel12 mon_zhat1(imon) = zhat1 mon_zhat2(imon) = zhat2 mon_phihat6(imon) = phihat6 mon_phihat12(imon) = phihat12 ; -- special for extra debugging END endif ; Determine the ML composition (comp=0 for group_num le 7, else comp=1) ; which is equivalent to comp=0 for mirror < 90, else comp=1 comp = 0; for Pt/C coating if mirror_arr(iphot) ge 90 then comp = 1; for W/Si coating angles = (180/!pi) * [ang1_arr(iphot),ang2_arr(iphot)]; a two element array ; converted from radians to degrees ;==================================================================================== ; ; What variations to include is determined by the variables ; use_as_coated, use_dmod_var, use_gamma_var, use_rough_var ; ; First the array 'params' is defined ; d_mod1 = 1.0 d_mod2 = 1.0 g_mod1 = 1.0 g_mod2 = 1.0 r_mod1 = 1.0 r_mod2 = 1.0 ; Locate the rows in the Master Table FINAL (MTF) w = where( mirror_arr(iphot) eq MTF_layer, nw ); mirror_arr(iphot) is current layer number uplo = MTF_uplo(w) u = where( uplo eq 'U' , nu ) l = where( uplo eq 'L' , nl ) ; There ought to be 12 or 24 of these (U + L) case nw of 12: begin ; use sextant iphi_rel = iphi_rel6 phihat = phihat6 mtf_index_u = w(u(sextant_num)) mtf_index_l = w(l(sextant_num)) ; -- special for extra debugging BEGIN if mon then mon_mtfu( imon ) = mtf_index_u ; -- special for extra debugging END end 24: begin ; use twelftant iphi_rel = iphi_rel12 phihat = phihat12 mtf_index_u = w(u(twelftant_num)) mtf_index_l = w(l(twelftant_num)) ; -- special for extra debugging BEGIN if mon then mon_mtfu( imon ) = mtf_index_u ; -- special for extra debugging END end else: begin print,'Bad nw' return end endcase ; Determine the recipe index (i.e. number - 1) recipe_idxu = MTF_recipe(mtf_index_u) - 1 recipe_idxl = MTF_recipe(mtf_index_l) - 1 ; Take action for very first layer that has recipe 0 (zero) ; ignore this for the time being if recipe_idxu lt 0 then recipe_idxu = 0 if recipe_idxl lt 0 then recipe_idxl = 0 ; Determine whether this is a DTU coated or RXO coated mirror rxo_coated_u = ( MTF_poschamb(mtf_index_u) eq 'rxo_pos' ) rxo_coated_l = ( MTF_poschamb(mtf_index_l) eq 'rxo_pos' ) idx_rxo = 0; -- choose a representative example if mon then mon_flag_rxo(imon) = 200*rxo_coated_u + rxo_coated_l mtfppu = MTF_pos_on_plate(mtf_index_u) mtfppl = MTF_pos_on_plate(mtf_index_l) mtfspu = MTF_separ(mtf_index_u) mtfspl = MTF_separ(mtf_index_l) ; ******* React if non-uniformity has been asked for ******************* if use_dmod_var or use_gamma_var or use_rough_var then begin ; IMark01 ; get modification of d-spacing based on NB's position dependent results ; but only if good results exist if use_dmod_var then begin ; IMark02 ; find d-spacing modifications for U section if rxo_coated_u then begin ; this is for an RXO mirror a = afit_rxo_t_coefs(*,idx_rxo) d_mod1 = ((a(0)*zhat1 + a(1))*zhat1 + a(2)) * ((a(3)*phihat + a(4))*phihat + a(5)) endif else begin ; this is for a DTU mirror w = where( mtfppu eq afit_t_pp and mtfspu eq afit_t_sp and comp eq afit_t_material, nw ) if nw eq 1 then begin a = afit_t_coefs(*,w(0)) d_mod1 = ((a(0)*zhat1 + a(1))*zhat1 + a(2)) * ((a(3)*phihat + a(4))*phihat + a(5)) endif else begin d_mod1 = 1.0 endelse endelse ; find d-spacing modifications for L section if rxo_coated_l then begin ; this is for an RXO mirror a = afit_rxo_t_coefs(*,idx_rxo) d_mod2 = ((a(0)*zhat2 + a(1))*zhat2 + a(2)) * ((a(3)*phihat + a(4))*phihat + a(5)) endif else begin ; this is for a DTU mirror w = where( mtfppl eq afit_t_pp and mtfspl eq afit_t_sp and comp eq afit_t_material, nw ) if nw eq 1 then begin a = afit_t_coefs(*,w(0)) d_mod2 = ((a(0)*zhat2 + a(1))*zhat2 + a(2)) * ((a(3)*phihat + a(4))*phihat + a(5)) endif else begin d_mod2 = 1.0 endelse endelse endif ; IMark02 if use_gamma_var then begin ; IMark03 ; find gamma modifications for U section but only if not RXO coated if not rxo_coated_u then begin w = where( mtfppu eq afit_g_pp and mtfspu eq afit_g_sp and comp eq afit_g_material, nw ) if nw eq 1 then begin a = afit_g_coefs(*,w(0)) g_mod1 = ((a(0)*zhat1 + a(1))*zhat1 + a(2)) * ((a(3)*phihat + a(4))*phihat + a(5)) endif else begin g_mod1 = 1.0 endelse endif ; find gamma modifications for L section but only if not RXO coated if not rxo_coated_l then begin w = where( mtfppl eq afit_g_pp and mtfspl eq afit_g_sp and comp eq afit_g_material, nw ) if nw eq 1 then begin a = afit_g_coefs(*,w(0)) g_mod2 = ((a(0)*zhat2 + a(1))*zhat2 + a(2)) * ((a(3)*phihat + a(4))*phihat + a(5)) endif else begin g_mod2 = 1.0 endelse endif endif ; IMark03 if use_rough_var then begin ; IMark04 ; find roughness modifications for U section but only if not RXO coated if not rxo_coated_u then begin w = where( mtfppu eq afit_r_pp and mtfspu eq afit_r_sp and comp eq afit_r_material, nw ) if nw eq 1 then begin a = afit_r_coefs(*,w(0)) r_mod1 = ((a(0)*zhat1 + a(1))*zhat1 + a(2)) * ((a(3)*phihat + a(4))*phihat + a(5)) endif else begin r_mod1 = 1.0 endelse endif ; find roughness modifications for L section but only if not RXO coated if not rxo_coated_l then begin w = where( mtfppl eq afit_r_pp and mtfspl eq afit_r_sp and comp eq afit_r_material, nw ) if nw eq 1 then begin a = afit_r_coefs(*,w(0)) r_mod2 = ((a(0)*zhat2 + a(1))*zhat2 + a(2)) * ((a(3)*phihat + a(4))*phihat + a(5)) endif else begin r_mod2 = 1.0 endelse endif endif ; IMark04 ; -- special for extra debugging BEGIN if mon then begin mon_dmod1( imon ) = d_mod1 mon_dmod2( imon ) = d_mod2 endif ; -- special for extra debugging END endif ; IMark01 ;====================== Using standard recipes ===================================================== if not use_as_coated then begin ; IMark05 ; get rcoefs from the standard recipe params1 = [recipe_std(recipe_idxu).N, recipe_std(recipe_idxu).dmin, $ recipe_std(recipe_idxu).dmax, recipe_std(recipe_idxu).c, $ recipe_std(recipe_idxu).topgam, recipe_std(recipe_idxu).gam] params1(1) = params1(1) * d_mod1 params1(2) = params1(2) * d_mod1 params1(4) = params1(4) * g_mod1 params1(5) = params1(5) * g_mod1 if rxo_coated_u then rough_factor = rrxorough else rough_factor = rsigma if bounce_arr(iphot) eq 3 or bounce_arr(iphot) eq 1 then begin coefu = get_rcoefs( earr, angles(0), rough*r_mod1*rough_factor, comp, params1 ) ; 'sigma' is entered as a scalar, same roughness throughout the stack endif else coefu = replicate(1.0, n_earr) if mon then mon_rco1(*,imon) = coefu(mon_idx) params2 = [recipe_std(recipe_idxl).N, recipe_std(recipe_idxl).dmin, $ recipe_std(recipe_idxl).dmax, recipe_std(recipe_idxl).c, $ recipe_std(recipe_idxl).topgam, recipe_std(recipe_idxl).gam] params2(1) = params2(1) * d_mod2 params2(2) = params2(2) * d_mod2 params2(4) = params2(4) * g_mod2 params2(5) = params2(5) * g_mod2 if rxo_coated_l then rough_factor = rrxorough else rough_factor = rsigma if bounce_arr(iphot) eq 3 or bounce_arr(iphot) eq 2 then begin coefl = get_rcoefs( earr, angles(1), rough*r_mod2*rough_factor, comp, params2 ) ; 'sigma' is entered as a scalar, same roughness throughout the stack endif else coefl = replicate(1.0, n_earr) if mon then mon_rco2(*,imon) = coefl(mon_idx) rcoef_sum = rcoef_sum + coefu*coefl endif ; IMark05 ;====================== Using as-coated recipes ===================================================== if use_as_coated then begin ; IMark06 ; get rcoefs from NB fittings to witness sample results ; but only if good results exist ; Locate the rows for the mirror shell for the current photon in the Master Table FINAL (MTF) w = where( mirror_arr(iphot) eq MTF_layer, nw ); mirror_arr(iphot) is current layer number witnesses = MTF_witness(w) uplo = MTF_uplo(w) u = where( uplo eq 'U' , nu ) l = where( uplo eq 'L' , nl ) ; There ought to be 12 or 24 of these (U + L) case nw of 12: begin ; use sextant witness_u = witnesses(u(sextant_num)); Name of witness for upper mirror witness_l = witnesses(l(sextant_num)); Name of witness for lower mirror if mon then mon_witness_u(imon) = witness_u if mon then mon_witness_l(imon) = witness_l if mon then mon_glass(imon) = sextant_num end 24: begin ; use twelftant witness_u = witnesses(u(twelftant_num)); Name of witness for upper mirror witness_l = witnesses(l(twelftant_num)); Name of witness for lower mirror if mon then mon_witness_u(imon) = witness_u if mon then mon_witness_l(imon) = witness_l if mon then mon_glass(imon) = twelftant_num end else: begin print,'Bad nw' return end endcase ; is there any match with Nicolai's results for the U mirror? wu = where( witness_u eq nb_sample, nwu ) if nwu eq 1 then begin ; yes, we have a match idxu = wu(0); Index in table of witness data params1 = [n(idxu), dmin(idxu)*d_mod1, dmax(idxu)*d_mod1, c(idxu), topgam(idxu)*g_mod1, $ gam(idxu)*g_mod1, gtot2(idxu)*g_mod1, gtot3(idxu)*g_mod1, gtot4(idxu)*g_mod1, $ gtot5(idxu)*g_mod1, $ dtot1(idxu), dtot2(idxu), dtot3(idxu), dtot4(idxu), dtot5(idxu), $ dtot6(idxu)] endif else begin ; calculate the standard coef idxu = -1 params1 = [recipe_std(recipe_idxu).N, recipe_std(recipe_idxu).dmin*d_mod1, $ recipe_std(recipe_idxu).dmax*d_mod1, recipe_std(recipe_idxu).c, $ recipe_std(recipe_idxu).topgam*g_mod1, recipe_std(recipe_idxu).gam*g_mod1] endelse if mon then if idxu ge 0 then mon_flag_usedfit(imon) = 200 ; is there any match with Nicolai's results for the L mirror? wl = where( witness_l eq nb_sample, nwl ) if nwl eq 1 then begin ; yes, we have a match idxl = wl(0) params2 = [n(idxl), dmin(idxl)*d_mod2, dmax(idxl)*d_mod2, c(idxl), topgam(idxl)*g_mod2, $ gam(idxl)*g_mod2, gtot2(idxl)*g_mod2, gtot3(idxl)*g_mod2, gtot4(idxl)*g_mod2, $ gtot5(idxl)*g_mod2, $ dtot1(idxl), dtot2(idxl), dtot3(idxl), dtot4(idxl), dtot5(idxl), $ dtot6(idxl)] endif else begin ; calculate the standard coef idxl = -1 params2 = [recipe_std(recipe_idxl).N, recipe_std(recipe_idxl).dmin*d_mod2, $ recipe_std(recipe_idxl).dmax*d_mod2, recipe_std(recipe_idxl).c, $ recipe_std(recipe_idxl).topgam*g_mod2, recipe_std(recipe_idxl).gam*g_mod2] endelse if mon then if idxl ge 0 then mon_flag_usedfit(imon) += 1 if bounce_arr(iphot) eq 3 or bounce_arr(iphot) eq 1 then begin if rxo_coated_u then rough_factor = rrxorough else rough_factor = rsigma if idxu lt 0 then begin coefu = get_rcoefs( earr, angles(0), rough*r_mod1*rough_factor, comp, params1 ) ; ---- Test for NaN returned wnan = where( coefu ne coefu, nwnan ) wf = where( 1 - finite(coefu), nwf ) if nwnan gt 0 or nwf gt 0 then begin print,'##1## NaN from get_rcoefs' print,format='("nwnan = ",i6," nwf = ",i6)', nwnan, nwf endif ; 'sigma' is entered as a scalar, same roughness throughout the stack endif else begin coefu = get_rcoefs( earr, angles(0), $ r_mod1*rough_factor*[rough_hz(idxu), rough_lz(idxu), 3.9, roughgrad(idxu)], $ comp, params1, cathodemod(*,idxu), interface=interface(idxu) ) ; 'sigma' is given as a 4-element array i.e. using layer individual roughnesses ; ---- Test for NaN returned wnan = where( coefu ne coefu, nwnan ) wf = where( 1 - finite(coefu), nwf ) if nwnan gt 0 or nwf gt 0 then begin print,'##2## NaN from get_rcoefs' print,format='("nwnan = ",i6," nwf = ",i6)', nwnan, nwf snapshot,v1=angles(0), v2=r_mod1*rough_factor*[rough_hz(idxu), rough_lz(idxu), 3.9, roughgrad(idxu)], $ v3=params1, t1='Angles', u1='deg', t2='Sigma', u2='A', t3='Params1', $ v4=r_mod1, t4='r_mod1' endif if mon then begin mon_flag_catdrop(imon) = 200*flag_catdrop(idxu) mon_flag_baddata(imon) = 200*flag_baddata(idxu) mon_flag_badfit(imon) = 200*flag_badfit(idxu) endif endelse endif else coefu = replicate(1.0, n_earr) if mon then mon_rco1(*,imon) = coefu(mon_idx) if bounce_arr(iphot) eq 3 or bounce_arr(iphot) eq 2 then begin if rxo_coated_l then rough_factor = rrxorough else rough_factor = rsigma if idxl lt 0 then begin coefl = get_rcoefs( earr, angles(1), rough*rough_factor*r_mod2, comp, params2 ) ; 'sigma' is entered as a scalar, same roughness throughout the stack ; ---- Test for NaN returned wnan = where( coefl ne coefl, nwnan ) wf = where( 1 - finite(coefl), nwf ) if nwnan gt 0 or nwf gt 0 then begin print,'##3## NaN from get_rcoefs' print,format='("nwnan = ",i6," nwf = ",i6)', nwnan, nwf endif endif else begin coefl = get_rcoefs( earr, angles(1), $ r_mod2*rough_factor*[rough_hz(idxl), rough_lz(idxl), 3.9, roughgrad(idxl)], $ comp, params2, cathodemod(*,idxl), interface=interface(idxl) ) ; 'sigma' is given as a 4-element array i.e. using layer individual roughnesses ; ---- Test for NaN returned wnan = where( coefl ne coefl, nwnan ) wf = where( 1 - finite(coefl), nwf ) if nwnan gt 0 or nwf gt 0 then begin print,'##4## NaN from get_rcoefs' print,format='("nwnan = ",i6," nwf = ",i6)', nwnan, nwf endif if mon then begin mon_flag_catdrop(imon) = mon_flag_catdrop(imon) + flag_catdrop(idxl) mon_flag_baddata(imon) = mon_flag_baddata(imon) + flag_baddata(idxl) mon_flag_badfit(imon) = mon_flag_badfit(imon) + flag_badfit(idxl) endif endelse endif else coefl = replicate(1.0, n_earr) if mon then mon_rco2(*,imon) = coefl(mon_idx) rcoef_sum = rcoef_sum + coefu*coefl endif ; IMark06 endfor; End loop over photons (rows in photon file) ;---------------------------------------------------------------------- ; End looping over all photons ---------------------------------- ;---------------------------------------------------------------------- fxaddpar,hdr,'DATE',ndate(3),'Date/time of file creation' fxaddpar,hdr,'TELESCOP','NuSTAR','Name of telescope (mission)' fxaddpar,hdr,'INSTRUME', mstr,'Name of instrument' fxaddpar,hdr,'FILE_IN',photon_file,'Name of input photon file' fxaddpar,hdr,'NROWS',nrows,'Number of photons in input file' fxaddpar,hdr,'MODE',mode,'Value of mode argument' fxaddpar,hdr,'RDMIN',rdmin,'dmin modification factor applied' fxaddpar,hdr,'RDMAX',rdmax,'dmax modification factor applied' fxaddpar,hdr,'RSIGMA',rsigma,'Roughness modification factor applied' fxaddpar,hdr,'RRXOROUG',rrxorough,'Roughness modification factor for RXO applied' fxaddpar,hdr,'PTOPTCST', getenv('PT_OPTCSTFILE'),'Opt. cst. used for Pt' fxaddpar,hdr,'COPTCST', getenv('C_OPTCSTFILE'),'Opt. cst. used for C' fxaddpar,hdr,'WOPTCST', getenv('W_OPTCSTFILE'),'Opt. cst. used for W' fxaddpar,hdr,'SIOPTCST', getenv('SI_OPTCSTFILE'),'Opt. cst. used for Si' ;;+ fxaddpar,hdr,'ROUGHNES',rough,'[angstrom] Value of roughness (sigma) parameter' fxaddpar,hdr,'ATTENUAT',atten,'Photon attenuation used' fxaddpar,hdr,'ATTENI',satten,'Photon attenuation subselection used' if keyword_set(bsel) then begin fxaddpar,hdr,'BSEL',bsel,'Bounce selection value' endif fxaddpar,hdr,'APERAREA',aperture_area,'[cm2] Aperture area' fxaddpar,hdr,'N_INJECT',n_inject,'Number of injected photons' fxaddpar,hdr,'N_TERMS',n_terms,'Number of photons included' fxaddpar,hdr,'SOFTWARE','calc_eff_area_nustar_v'+version+'.pro','IDL generating program' fxaddpar,hdr,'RESPONSI','Niels J. Westergaard','Responsible person' fxaddpar,hdr,'TUNIT1','keV','Unit of first column (ENERGY)' fxaddpar,hdr,'TUNIT2','cm2','Unit of second column (EFF_AREA)' fxaddpar,hdr,'COMMENT','The effective area for the optic (FM1 or FM2) as calculated' fxaddpar,hdr,'COMMENT','using Dave Windts IMD software. For mode xxx1 the measured' fxaddpar,hdr,'COMMENT','witness sample fitted coating characteristics as of 2012-01-20' fxaddpar,hdr,'COMMENT','by Nicolai Brejnholt have been used.' fxaddpar,hdr,'COMMENT','mode = 0001: Using as-coated (in contrast to nominal recipes)' fxaddpar,hdr,'COMMENT','mode = 0010: Varying d-spacing' fxaddpar,hdr,'COMMENT','mode = 0100: Varying gamma values' fxaddpar,hdr,'COMMENT','mode = 1000: Varying roughness values' wrmfitscols,outname,"ENERGY",earr,"EFF_AREA",atten*rcoef_sum*aperture_area/n_inject, $ "AVG_RCOEF",rcoef_sum/n_terms,kwds=hdr print,'Has written: '+outname ; Only output the monitor file for as-coated cases if mon and use_as_coated then begin for ie = 0, 6 do begin fxaddpar,hdr,'ENERGY0'+stri(ie+1),earr(mon_idx(ie)),'[keV] Energy for RCOEFs' endfor wrmfitscols,filemon, "LAYER", mon_layer, "GLASSPIECE", mon_glass, "WITNESS_U", $ mon_witness_u,"WITNESS_L", mon_witness_l, "SECTOR", mon_s_num, "SUBSECTOR", mon_t_num, $ "I1Z", mon_i1z, "I2Z", mon_i2z, "Z1_REL", mon_z1_rel, $ "ANGLE_IN1", mon_ang1, "ANGLE_IN2", mon_ang2, "RCOEF1", mon_rco1, "RCOEF2", mon_rco2, $ "PHI", mon_phi, $ "PHI_REL6", mon_phi_rel6, "PHI_REL12", mon_phi_rel12, $ "MTFU_IDX", mon_mtfu, "DMOD1", mon_dmod1, "DMOD2", mon_dmod2, $ "ZHAT1", mon_zhat1, "ZHAT2", mon_zhat2, "PHIHAT6", mon_phihat6, "PHIHAT12", mon_phihat12, $ "FLAG_CATDROP", mon_flag_catdrop, "FLAG_BADDATA", mon_flag_baddata, $ "FLAG_BADFIT", mon_flag_badfit, "FLAG_RXO", mon_flag_rxo, "FLAG_USEDFIT", mon_flag_usedfit, kwds=hdr print,'Has written: '+filemon endif else print,'Skipped writing monitor file ...' end