pro for_uv_densminmax,ObsPramsStruct,ModPramsStruct,rpos,chromoheight,theta3d_chromo=theta3d_chromo,phi3d_chromo=phi3d_chromo,chromopatch=chromopatch,do_fullsphere=do_fullsphere,fullsphere_dmin=fullsphere_dmin,fullsphere_dmax=fullsphere_dmax,nreinit=nreinit,nowidgmess=nowidgmess,ARw=ARw ;------------------------------------------------------------------------ ; Purpose: ; ; This routine figures out the min and max density (or ion density) ; to be used in weighting PSI active regions for anisotropic illumination ; from below a scattering point -- and uses to set AR weight ; *special cse chromopatch = 0 will instead replace input ; *keyword chromorad with self-consistent collisional intensity ; *evaluated at height chromoheight ; ; Inputs: ; ; ObsPramsStruct, ModPramsStruct ; ; RPOS - needed for FOR_IONDENS input keyword ; ; CHROMOHEIGHT -- height of base of scattering cone where the source emission is located ; **for now assuming single number ; ; KEYWORD inputs ; ; THETA3D_CHROMO, PHI3D_CHROMO 3D global coordinates of cone footpoints ; not used first time called by FOR_INTENSINT global calculation ; ; CHROMOPATCH -- see FOR_OBSDEFAULTS for explanation ; but note if it set less than zero, this subroutine wont be called (instead FOR_ARWEIGHT) ; **If it is zero, ; **the ion contribution to intensity will be evaluated at the footpoints ; **and the passed back quantity will be dimensionless ARw=I_footpoints/chromorad ; **since it will cancel chromorad dependency in for_uvmodel (within ionemiss term) ; ; DO_FULLSPHERE is _not_ a global FORWARD keyword ; evaluates minmax below with a full sphere THETA, PHI array ; intead of THETA3D/PHI3D_CHROMO ; --> output FULLSPHERE_DMIN/DMAX ; =1 set in FOR_INTENSINT ; and used to do an initial run (before doing all the cone calculations ; for scattering points) to establish min/max ion density to scale ; the chromorad at cone footpoints** ; =0 when the calculation is being done for cone footpoints ; as called by FOR_UV_STOKES ; ; FULLSPHERE_DMIN/DMAX also not a global FORWARD keyword ; used to scale the CHROMOPATCH with a consistent global scaling ; If DO_FULLSPHERE = 1 (i.e., the first call from FOR_INTENSINT) ; these are the output of this code (DMIN/DMAX) ; will be used in the later calls from FOR_UV_STOKES ; If DO_FULLSPHERE = 0 (i.e., these later calls from FOR_UV_STOKES ; now these are inputs, ; will be used to scale the footpoint iondensity weighting ; **UNLESS! CHROMOPATCH=0 -- in which case, ; *CHROMORAD will be self-consistently replaced by intensity ; *see below) ; ; NREINIT -- needed for call to numerical models ; NOWIDGMESS -- needed for call to subroutines ; ; KEYWORD Outputs: ; ; ARw -- the array of weights ; equal to 1 everywhere except points above 500000 threshold ; and there it is set to values ranging from 1 to 1 + chromopatch ; **(unless chromopatch = 0, see above) ; note -- ignored and not passed back if DO_FULLSPHERE is set (=1) ; ; FULLSPHERE_DMIN/DMAX -- ; outputs if the code is called with DO_FULLSPHERE=1 (e.g.. by FOR_INTENSINT) ; ; Called by FOR_UV_STOKES, FOR_INTENSINT ; Calls PSIMAS, FOR_IONDENS ;------------------------------------------------------------------------ ; ; July 2025 -- added dowidgmess keyword and updates ; to enable other codes besides PSIMAS ; and DO_FULLSPHERE=-1 option to forward model ; intensity at footpoints, not just scaling to ion density ; also added Lya as an option in call to for_iondens ; and removed conditional for "other lines" ; should be able to do them all, even if not part of original EIS set ; also removed the calculation of ObsPramsStructChromo-- instead ; now have edited for_idset so it can handle that the UV cases ; have line=i,q,u,v ; August 2025/February 2026 -- updated and cleaned up comments ; March 2026 -- vectorized grid setup, COMMON-block collision rate cache, ; CASE instrument dispatch, unique-temperature interpolation ; RK - COMMON block caches collision rate table across calls so that the ; rate file is read from disk only once per instrument, not once per ; LOS cone. cached_instrument guards against instrument changes ; within a single FORWARD run. COMMON uv_collrates_cache, cached_instrument, cached_t, cached_c usewindows=0 if strupcase(!version.os_family) eq 'WINDOWS' then usewindows=1 if keyword_set(do_fullsphere) eq 1 then begin ; Makes theta3d_chromo, phi3d_chromo with global grid. ; This should only be set when called from for_intensint. ; ; RK - replaced two scalar for-loops with fully vectorized rebin calls. theta3d_chromo = rebin(reform(dindgen(180)*!dtor, 1, 180), 360, 180) phi3d_chromo = rebin(dindgen(360)*!dtor, 360, 180) print,'**************************' print,'**************************' print,'SETTING UP FULL SPHERE ION DENSITY MIN MAX' print,'**************************' print,'**************************' endif else begin ; do_fullsphere=0 - if chromopatch ne 0 then begin ; calculate weights ; print,'**************************' ; print,'USING FULL SPHERE ION DENSITY MIN MAX FOR CHROMORAD WEIGHTS' ; print,'**************************' endif else begin ; print,'**************************' ; print,'**************************' ; print,'CALCULATING MODEL CHROMORAD INTENSITY $ ; AT FOOTPOINTS FOR SELF-CONSISTENT CHROMORAD WEIGHTS' ; print,'**************************' ; print,'**************************' endelse endelse ; RK - replicate() allocates and fills directly instead of ; constructing an intermediate zero-valued double array via *0.d0 chromoheight_2d = replicate(double(chromoheight), size(theta3d_chromo, /dimensions)) ; Note CHROMOHEIGHT LT 1.01 IS UNLIKELY TO DO ANY WEIGHTING for PSI ; because of condition placed on temperature to remove overbroad transition region IF strupcase(ModPramsStruct.Name) EQ 'NUMCUBE' $ or strupcase(ModPramsStruct.Name) EQ 'ADAPTCUBE' $ or strupcase(ModPramsStruct.Name) EQ 'AWSOM' $ or strupcase(ModPramsStruct.Name) eq 'PSIMAS' THEN BEGIN IF strupcase(ModPramsStruct.Name) EQ 'NUMCUBE' THEN $ call_procedure,ModPramsStruct.Name,chromoheight_2d,theta3d_chromo,phi3d_chromo,ModPramsStruct,ModSolStructChromo,nreinit=nreinit,nowidgmess=nowidgmess $ else call_procedure,ModPramsStruct.Name,chromoheight_2d,theta3d_chromo,phi3d_chromo,ModPramsStruct,ModSolStructChromo,nreinit=nreinit endif else call_procedure,ModPramsStruct.Name,chromoheight_2d,theta3d_chromo,phi3d_chromo,ModPramsStruct,ModSolStructChromo ; ModSolStructChromo is the output of this procedure for_iondens,ObsPramsStruct,ModSolStructChromo,ModPramsStruct,chromoheight,rpos,iondens,nowidgmess=nowidgmess,nochromoadd=1 ; NOTE NOCHROMOADD IS NOT A GLOBAL KEYWORD ; forced here so that we dont add the background to the model iondens ; Now calculate collision rates ; collisional rates (from CHIANTI) ; RK - instrument string is uppercased once and reused throughout. ; CASE replaces a chain of sequential IF statements ; raises an explicit error for unknown instruments ; instead of silently leaving collfile undefined. inst = strupcase(ObsPramsStruct.instrument) colldir = concat_dir(GET_ENVIRON('FORWARD'), 'DEFAULTS') CASE inst OF 'OVI1032': collfile = concat_dir(colldir, 'OVI1032_rates.txt') 'OVI1037': collfile = concat_dir(colldir, 'OVI1037_rates.txt') 'LYA': collfile = concat_dir(colldir, 'LYA_rates.txt') 'NEVIII770': collfile = concat_dir(colldir, 'NEVIII770_rates.txt') 'NEVIII780': collfile = concat_dir(colldir, 'NEVIII780_rates.txt') 'MGIX706': collfile = concat_dir(colldir, 'MGIX706_rates.txt') ELSE: message, 'for_uv_densminmax: unknown instrument for collision rates: ' + inst ENDCASE if usewindows eq 1 then collfile = str_replace(collfile, '/', '\') ; RK - collision rate table is read from disk only on the first call, ; or when the instrument changes. All subsequent calls for the same ; instrument retrieve t and c directly from the COMMON block, eliminating ; repeated file I/O across every LOS cone call from FOR_UV_STOKES. if n_elements(cached_instrument) eq 0 then begin readcol, collfile, t, c, /silent cached_t = t cached_c = c cached_instrument = inst endif else if cached_instrument ne inst then begin readcol, collfile, t, c, /silent cached_t = t cached_c = c cached_instrument = inst endif else begin t = cached_t c = cached_c endelse ; RK - interpolate onto unique log-temperatures only, then scatter ; back to the full grid using value_locate. For smooth MHD/PFSS models ; where large regions share nearly identical temperatures this avoids ; redundant interpol evaluations across the full 2-D array. logT = alog10(ModSolStructChromo.Temp) logT_flat = logT[*] sort_idx = sort(logT_flat) uniq_idx = uniq(logT_flat, sort_idx) uniq_logT = logT_flat[uniq_idx] uniq_crates = interpol(c, t, uniq_logT) ; value_locate returns the index of the nearest lower bound ; in the sorted unique array for every element of logT_flat crates = reform(uniq_crates[value_locate(uniq_logT, logT_flat)], $ size(logT, /dimensions)) ; units of cm^3 s^-1 RSun_cm = 6.95700d+10 ;solar radius in cm testtrans=where(finite(iondens) ne 0 and iondens ne 0.d0 $ and finite(ModSolStructChromo.Dens) ne 0.d0 and $ ModSolStructChromo.Dens ne 0.d0, ntt) ; set min/max to global values ; will not exist first time called from FOR_INTENSINT ; and will be created as input ; but then will for all the calls from FOR_UV_STOKES ; so will be used as input if exist(fullsphere_dmin) then mindens=fullsphere_dmin else $ mindens=min(iondens[testtrans]) if exist(fullsphere_dmax) then maxdens=fullsphere_dmax else $ maxdens=max(iondens[testtrans]) if do_fullsphere eq 0 then begin ; ion density at footpoints for calls from FOR_UV_STOKES I_footpoints=iondens*0.d0 densuse=iondens[testtrans] if ntt ne 0 and maxdens ne mindens then begin ; ; FIRST possibility: chromopatch gt 0 ; will uses *observed* chromorad value (applied in for_uv_model.pro) ; (*observed = e.g. from SUMER, set in for_specdefaults) ; and applies chromopatch as weight to ion density at the emitting height chromoheight if chromopatch ne 0 then begin ARw[testtrans] = 1.d0 + chromopatch*(densuse - mindens)/(maxdens - mindens) endif else begin ; ; SECOND possibility: chromopatch = 0 ; substitute a self-consistent collisional intensity at the footpoints ; that will cancel the *observed* chromorad that is introduced in for_uvmodel ; (part of the ionemiss term that multipies the Stokes vector) ; ; chromorad has units of photons cm^-2 s^-1 sr^-1 ; because it is the total intensity (integrated over solar disk) ; ; see for_uvmodel.pro for collision calculation ; units of ergs cm^-3 s^-1 sr^-1 ; so collemiss/photon_energy*RSun_cm will have units photons cm^2 s^-1 sr^-1 ; I_footpoints[testtrans]=0.83*RSun_cm*densuse*ModSolStructChromo.Dens[testtrans]*crates[testtrans]/(4.*!pi) ; need to do something with the points that intersect T < threshold ; here we set to the minimum intensity for the rest of the cone footpoints (satisfying testtrans) ; might make more sense to set to zero ; ; minIfoot=min(I_footpoints[testtrans]) minIfoot=0.d0 ARw=minIfoot/ObsPramsStruct.SpecPrams.ChromoRad ARw[testtrans] = I_footpoints[testtrans]/ObsPramsStruct.SpecPrams.ChromoRad endelse test = where (finite(arw) eq 0) if min(test) ne -1 then stop endif else begin ; ; RK - No valid TR footpoints (ntt=0) or all footpoints have identical density (mindens=maxdens). ; For chromopatch ne 0: ARw stays at 1 (uniform observed chromorad -- correct default). ; For chromopatch eq 0: ARw must be zeroed. The self-consistent model has no collisional ; emission at these footpoints, so no radiation illuminates the scattering cone from below. ; Leaving ARw=1 here would incorrectly apply the full observed chromorad (e.g. SUMER value) ; to directions where the model predicts no TR emission -- producing spurious bright features ; at coronal hole footpoints. ; SG - Agreed. This is consistent with the decision above for points T < threshold ; -- we considered setting to minimum intensity of rest of cone but ended up setting to zero ; if chromopatch eq 0 then ARw = ARw * 0.d0 endelse endif else begin ; set these as outputs fullsphere_dmin=mindens fullsphere_dmax=maxdens print,'fullsphere_dmin=',mindens print,'fullsphere_dmax=',maxdens endelse end