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 -- height of base of scattering cone where the source emission is located ; single number, default r=1 ; (may vary with theta,phi in future, but this should *NOT* be done vs ; scattering point theta, phi, rather cone footpoints, which will be a bit trickier) ; ; CHROMO_THMIN, CHROMO_THMAX, CHROMO_PHMIN, CHROMO_PHMAX -- ; patch of anisotropic (active region) radiation ; Default unset -- all set to 0 ; CHROMOPATCH -- value to multiply by, default 10 ; 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 ;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) 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) if ntst ne 0 then stop ;------------------------------------------------------------------------- ; 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 ; PATCH: ; IF WANT TO COMPARE CHROMO ANISOTROPIC CASE TO CASE NOT ANISOTROPIC - ; first SET CHROMOPATCH=10 (or whatever) and then run again with CHROMOPATCH=0 if strupcase(ModPramsStruct.Name) ne 'PSIMAS' then ARw = for_arweight(chromoheight,theta3d_chromo,phi3d_chromo,chromo_thmin,chromo_thmax,chromo_phmin,chromo_phmax,chromopatch) else begin ARw = theta3d_chromo*0.d0 + 1.d0 ; PSI case if chromoheight gt 1.d0 and chromopatch ne 0.d0 then begin for_uv_densminmax,ObsPramsStruct,ModPramsStruct,rpos,chromoheight,theta3d_chromo=theta3d_chromo,phi3d_chromo=phi3d_chromo,chromopatch=chromopatch,ARw=Arw,fullsphere_dmin=fullsphere_dmin,fullsphere_dmax=fullsphere_dmax,nreinit=nreinit 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 end