pro for_uv_stokes,ObsPramsStruct,ModPramsStruct,r3D,theta3D,phi3D,thmod,phmod,thetaB,chiB,bigtheta,bigchi,Bmag,Vwind,w_par,aniso,nangleint,w_d,blend,einsteinA,gJ,chromoheight,chromo_thmin,chromo_thmax,chromo_phmin,chromo_phmax,chromopatch,ulimb,instrument,StokesUV_I,StokesUV_Q,StokesUV_U,nreinit=nreinit,nowidgmess=nowidgmess,fullsphere_dmin=fullsphere_dmin,fullsphere_dmax=fullsphere_dmax ;------------------------------------------------------------------------ ; ; Modified from pol_dim1.pro ; Silvano Fineschi ; for_uv_stokes.pro ;------------------------------------------------------------------------ ; Purpose: ; ; This routine computes the integrals of the Stokes parameters over the ; azimuthal and polar angles of the propagation vector of the chromospheric ; radiation coming from the spherical cap of the chromospere limited by a ; cone centered on the scattering point. ; ; Called by FOR_UVMODEL ; Calls FOR_D, FOR_P2_00, FOR_P2_10, FOR_P2_20, FOR_GLOBAL_LOCAL_ROTATE, FOR_ARWEIGHT,FOR_IONDENS,FOR_ARWEIGHT, PSIMAS,FOR_UV_DENSMINMAX ;------------------------------------------------------------------------ ; ; Inputs: ; All are points in this version ; ; R3D, radial (z') heliocentric height [Rsun] ; THETA3D,PHI3D spherical angular coordinates of scatterer position - radians ; **OBSERVER FRAME OF REFERENCE** ; --needed for calculating scattering angles ; THMOD, PHMOD spherical coordinates of scatterer position - radians ; **MODEL FRAME OF REFERENCE** ; *(rotated by B angle and central meridian and potentially thetao)* ; --needed for calculating non-isotropic radiation from below (ARWeight) ; ; (B vector spherical coordinated in the local solar vertical frame) ; THETAB, CHIB: polar and azimuthal angles of magnetic vector [rad] ; BIGTHETA: polar angle between local solar vertical and the line of sight [rad] ; BIGCHI: angle between the plane containing the ; line of sight and plane perpendicular to the plane of the sky ; BMAG: magnitude of magnetic vector [gauss] ; ; VWIND -- magnitude of solar wind speed (direction assumed to be along field) [km/s] ; W_PAR -- 1/e Gaussian width for velocity distribution of the coronal ion parallel the magnetic field [km/s] ; ANISO -- parameter describing anisotropy of coronal velocity distribution: ; 1/e Gaussian width for velocity distribution of the coronal ion perpendicular to the magnetic field ; W_PERP=W_PAR*ANISO; [1 = isotropic distribution] ; NANGLEINT -- how many steps to take in solid angle integration ; ; W_D -- 1/e Gaussian width for velocity distribution of the transition region ion [km/s] ; BLEND - line blending factor due to transition ; EINSTEINA: Einstein coefficient for spontaneous emission [1/s] ; GJ: Landè factor for the atomic transition ; ; CHROMOHEIGHT, CHROMO_THMIN, CHROMO_THMAX, CHROMO_PHMIN, CHROMO_PHMAX,CHROMOPATCH ; see explanation in FOR_OBSDEFAULTS ; ; ULIMB -- limb darkening ; INSTRUMENT ; ; Keyword Inputs ; ; FULLSPHERE_DMIN, FULLSPHERE_DMAX -- if set, will use in lieu of ; getting density min/max for cone footpoints if PSI ; ; NREINIT - to avoid reloading the PSI cube into memory ; NOWIDGMESS - to suppress popup widgets ; Outputs: ; ; STOKESUV_I,STOKESUV_Q,STOKESUV_U ; ; Nov 2024 -- added hooks for limb brightening functions -- RK ; Dec 2024 -- passed through chromoheight ; Spring 2025 -- added AR weight from below and ability to run psi for below ; June 2025 -- ; Changed asin/acos to ratio of floats to avoid numerical error when ; numerical nearly = denominator ; updated for fullsphere_dmin/dmax ; Aug 2025 -- commented out check for chromoheight gt 1 because not sure why there ; also put in fix for scatter below emitter ; Jan 2026 -- cleaned up comments etc on CHROMO* ; also fixed bug where CHROMOPATCH < 1 ; resulted in negative intensity ;Hanle effect parameter => gamma = 0.88 * B * g(J')/A(J',J) ; where B is in gauss and A(J',J) is 1e.7 s units ; See Eq. 16 in Landi Degl'Innocenti, Solar Phys. 102 (1-20) (1985) gamma = 0.88 * Bmag * gJ/(einsteinA/1.e7) ; Geometrical parameters angles={theta:bigtheta,chi:bigchi,thetaB:thetaB,chiB:chiB} ; Doppler dimming parameters widths={vwind:vwind,w_d:w_d,w_perp:aniso*w_par,w_par:w_par} ; polar angle subtending the solar disk ;theta_limb=asin(chromoheight/r3D) ; *note if we change chromoheight to vary with cone footpoints ; *the issue will be if it has a different dimension from r3D ; * which has the dimension of the scatteirng points, not the footpoints theta_limb=asin(float(chromoheight)/float(r3D)) ; ; THETA_LIMB is ; polar angle limit of the solid angle subtending the Sun from the heliocentric height "R3D" [rad] ; ; this next was for diagnosing weird numerical error where ; identical chromoheight, r3d threw a NaN ; now using FLOAT above to avoid this test = where(finite(theta_limb) ne 1,ntst) ; temporary fix for situation where chromoheight gt scattering point if ntst ne 0 then begin print,'scatter point below emitter; moving emitter down to scatterer' theta_limb=1.d0 ;stop endif ;------------------------------------------------------------------------- ; Main program ;------------------------------------------------------------------------- ; Integration over the spherical coordinates of the incoming radiation ; The discrete angular steps dtheta and dchi are defined for a reasonable integration over the spherical coordinates of the incoming radiation nsteps=nangleint dtheta_prime=theta_limb/nsteps dchi_prime=2.*!dpi/nsteps ;creating the angular vectors theta_prime=dindgen(nsteps)*dtheta_prime chi_prime=dindgen(nsteps)*dchi_prime ; Converting the angular vectors in 2D matrices in the 2D space of dtheta and dchi theta_prime_2D=theta_prime # replicate(1.,nsteps) chi_prime_2D=replicate(1.,nsteps) # chi_prime ; mu_2D=cos(zeta), where zeta is angle between solar vertical and the ray towards the scatterer at the point where it exits the surface (Casini and Judge 1999 Fig. 4) rpos= sin(theta_prime_2D)/sin(theta_limb) mu_2D = sqrt(1.-(rpos)^2) allen_limb_2D = mu_2D * 0. + 1. ; 2D Center-to-limb variation from Allen function if strpos(strupcase(instrument),'NEVIII') ne -1 and ulimb ne 0 then $ allen_limb_2D = ulimb * allen_fun_ne8(mu_2D) if strpos(strupcase(instrument),'OVI') ne -1 and ulimb ne 0 then $ allen_limb_2D = ulimb * allen_fun_o6(mu_2D) ; 2D Doppler dimming matrix D_2D = for_d(theta_prime_2D, chi_prime_2D, angles, widths) * sin(theta_prime_2D)*dtheta_prime*dchi_prime/(4.*!pi) ; ; figure out coordinates to compare to nonisotropic radiation from below ; for_global_local_rotate,chromoheight,theta_prime_2D,chi_prime_2D,r3D,thmod,phmod,$ theta3d_chromo,phi3d_chromo ; establish weights that will act as multiplier on NORMCHROMORAD in FOR_UVMODEL ; CHROMOPATCH: ; if set gt 0 will calculate weights based on iondens ; if set eq 0 will calculate self-consistent weights based on intensity ; if set lt 0 then will make a patch multiplier= abs(chromopatch) ; chromopatch=-1 will not produce a patch (multiply by 1) ; and will only do limb darkening if chromopatch lt 0 then ARw = for_arweight(chromoheight,theta3d_chromo,phi3d_chromo,chromo_thmin,chromo_thmax,chromo_phmin,chromo_phmax,abs(chromopatch)) else begin ARw = theta3d_chromo*0.d0 + 1.d0 ; if chromoheight gt 1.d0 then begin for_uv_densminmax,ObsPramsStruct,ModPramsStruct,rpos,chromoheight,theta3d_chromo=theta3d_chromo,phi3d_chromo=phi3d_chromo,chromopatch=chromopatch,ARw=ARw,do_fullsphere=0,fullsphere_dmin=fullsphere_dmin,fullsphere_dmax=fullsphere_dmax,nreinit=nreinit,nowidgmess=nowidgmess ; endif endelse StokesUV_I_2D=(1.+blend*for_p2_00(theta_prime_2D, chi_prime_2d, angles, gamma)) * D_2D * allen_limb_2D * ARw StokesUV_Q_2D=blend*for_p2_10(theta_prime_2D, chi_prime_2d, angles, gamma) * D_2D * allen_limb_2D * ARw StokesUV_U_2D=blend*for_p2_20(theta_prime_2D, chi_prime_2d, angles, gamma) * D_2D * allen_limb_2D * ARw ;------------------------------------------------------------------------ ; Results Output ;------------------------------------------------------------------------ StokesUV_I=total(StokesUV_I_2D) StokesUV_Q=total(StokesUV_Q_2D) StokesUV_U=total(StokesUV_U_2D) if finite(StokesUV_I) eq 0. then begin stop print,'StokesI',(StokesUV_I) print,'StokesI_2D',minmax(StokesUV_I_2D) print,'ARw',minmax(ARw) endif ;print,'ARw',minmax(ARw) ;print,'chromopatch',chromopatch ;print,'chromoheight',chromoheight end