subroutine integrate_emission_spherical(igrid,numXI1,numXI2,xI1,xI2,dxI,fl,&
datatype,EM)
integer, intent(in) :: igrid,numXI1,numXI2
double precision, intent(in) :: xI1(numXI1),xI2(numXI2)
double precision, intent(in) :: dxI
type(te_fluid), intent(in) :: fl
character(20), intent(in) :: datatype
double precision, intent(inout) :: EM(numXI1,numXI2)
integer :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,ixO1,ixO2,&
ixO3,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,ix1,ix2,ix3,i,j
double precision, allocatable :: flux(:,:,:),Ne(:,:,:)
integer :: ixPmin1,ixPmin2,ixPmin3,ixPmax1,ixPmax2,ixPmax3,ixP1,ixP2,&
ixP3,nSubC1,nSubC2,nSubC3,iSubC1,iSubC2,iSubC3
double precision :: xSubP1,xSubP2,dxSubP,xerfmin1,xerfmin2,xerfmin3,&
xerfmax1,xerfmax2,xerfmax3,fluxsubC,RsubC
double precision :: TBsubC,PBsubC
double precision :: xSubC(1:3),dxSubC1,dxSubC2,dxSubC3,xCent(1:2),&
xSubC_car(1:3)
double precision :: R_thick,dotp,dvolume,R_occult,Rc
double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
logical :: sun_back_side,emit
integer :: mass
double precision :: logTe
character (30) :: ion
double precision :: lineCent
double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
double precision :: RHESSI_rsl,area_1AU,arcsec,pixel
ixOmin1=ixmlo1;ixOmin2=ixmlo2;ixOmin3=ixmlo3;
ixOmax1=ixmhi1;ixOmax2=ixmhi2;ixOmax3=ixmhi3;
ixImin1=ixglo1;ixImin2=ixglo2;ixImin3=ixglo3;
ixImax1=ixghi1;ixImax2=ixghi2;ixImax3=ixghi3;
if (SI_unit) then
arcsec=7.25d5/unit_length
else
arcsec=7.25d7/unit_length
endif
allocate(flux(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3))
if (datatype=='image_euv') then
! get local EUV flux and velocity
call get_EUV(wavelength,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,ps(igrid)%w,&
ps(igrid)%x,fl,flux)
flux(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)=flux(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)/instrument_resolution_factor**2 !adjust flux due to artifical change of resolution
call get_line_info(wavelength,ion,mass,logTe,lineCent,spaceRsl,wlRsl,&
sigma_PSF,wslit)
pixel=spaceRsl*arcsec
sigma0=sigma_PSF*pixel
else if (datatype=='image_sxr') then
!get local SXR flux photons cm-3,cm-3,cm-3 s-1,s-1,s-1 (cgs) or photons m-3,m-3,m-3 s-1,s-1,s-1 (SI)
call get_SXR(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,ixOmin1,&
ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,ps(igrid)%w,ps(igrid)%x,fl,&
flux,emin_sxr,emax_sxr)
RHESSI_rsl=2.3d0/instrument_resolution_factor
sigma_PSF=1.d0
pixel=RHESSI_rsl*arcsec
sigma0=sigma_PSF*pixel
area_1AU=2.81d27
endif
! integrate emission
R_thick=R_opt_thick*const_Rsun/unit_length
do ix1=ixOmin1,ixOmax1
do ix2=ixOmin2,ixOmax2
do ix3=ixOmin3,ixOmax3
x_sph(1:3)=ps(igrid)%x(ix1,ix2,ix3,1:3)
dx_sph(1:3)=ps(igrid)%dx(ix1,ix2,ix3,1:3)
dxl(1)=dx_sph(1) ! cell size in length
dxl(2)=x_sph(1)*dx_sph(2) ! cell size in length
dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3) ! cell size in length
!dividing a cell to several sub-cells to get more accurate integrating values
nSubC1=1;nSubC2=1;nSubC3=1;
call get_unit_vector_spherical(x_sph,unitv_r,unitv_theta,unitv_phi)
call dot_product_loc(unitv_r,vec_xI1,dotp)
nSubC1=max(nSubC1,ceiling(dxl(1)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_r,vec_xI2,dotp)
nSubC1=max(nSubC1,ceiling(dxl(1)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_theta,vec_xI1,dotp)
nSubC2=max(nSubC2,ceiling(dxl(2)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_theta,vec_xI2,dotp)
nSubC2=max(nSubC2,ceiling(dxl(2)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_phi,vec_xI1,dotp)
nSubC3=max(nSubC3,ceiling(dxl(3)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_phi,vec_xI2,dotp)
nSubC3=max(nSubC3,ceiling(dxl(3)*abs(dotp)/(dxI/2.d0)))
! integrate sub-cells
do iSubc1=1,nSubC1
! sub-cell center coordinate in spherical
xSubC(1)=x_sph(1)-half*dx_sph(1)+(iSubC1-half)*dx_sph(1)/nSubC1
RsubC=xSubC(1)
dxSubC1=dx_sph(1)/nSubC1 ! sub-cell size in length
do iSubc2=1,nSubC2
! sub-cell center coordinate in spherical
xSubC(2)=x_sph(2)-half*dx_sph(2)+(iSubC2-half)*dx_sph(2)/nSubC2
dxSubC2=xSubC(1)*dx_sph(2)/nSubC2 ! sub-cell size in length
dxSubC3=xSubC(1)*dsin(xSubC(2))*dx_sph(3)/nSubC3 !sub-cell size in length
dvolume=dxSubC1*dxSubC2*dxSubC3
if (datatype=='image_euv') then
if (SI_unit) then
fluxSubC=flux(ix1,ix2,ix3)*dvolume*unit_length*1.d2/dxI/dxI !DN s-1,s-1,s-1
else
fluxSubC=flux(ix1,ix2,ix3)*dvolume*unit_length/dxI/dxI !DN s-1,s-1,s-1
endif
else if (datatype=='image_sxr') then
!sub-cell SXR flux at 1 AU [photons s-1,s-1,s-1 cm-2],cm-2],cm-2]
fluxSubC=flux(ix1,ix2,ix3)*dvolume*unit_length**3/area_1AU
endif
! enter integration if flux large enough
if (fluxSubC>smalldouble) then
do iSubc3=1,nSubC3
! sub-cell center coordinate in spherical
xSubC(3)=x_sph(3)-half*dx_sph(3)+&
(iSubC3-half)*dx_sph(3)/nSubC3
call get_cor_image_spherical(xSubC,xCent)
Rc=dsqrt(xCent(1)**2+xCent(2)**2) !distance to sun center (on the image plane)
!
! whether the local emitted photons can arrive the telescope
call spherical_to_cartesian(xSubC,xSubC_car)
call dot_product_loc(vec_LOS,xSubC_car,dotp)
sun_back_side=.true.
if (dotp<0.d0) sun_back_side=.false.
! whether the local emission can reach the telescope
if (sun_back_side) then
emit=.false.
if (Rc>R_thick) emit=.true.
else
emit=.true.
if (xSubC(1)<=R_thick) emit=.false.
endif
!
if (emit) then
! mapping the 3D coordinate to location at the image
! distribution at nearby pixels
ixP1=floor((xCent(1)-(xI1(1)-half*dxI))/dxI)+1
ixP2=floor((xCent(2)-(xI2(1)-half*dxI))/dxI)+1
ixPmin1=max(1,ixP1-3)
ixPmax1=min(ixP1+3,numXI1)
ixPmin2=max(1,ixP2-3)
ixPmax2=min(ixP2+3,numXI2)
do ixP1=ixPmin1,ixPmax1
do ixP2=ixPmin2,ixPmax2
xerfmin1=((xI1(ixP1)-half*dxI)-&
xCent(1))/(sqrt(2.d0)*sigma0)
xerfmax1=((xI1(ixP1)+&
half*dxI)-xCent(1))/(sqrt(2.d0)*sigma0)
xerfmin2=((xI2(ixP2)-half*dxI)-&
xCent(2))/(sqrt(2.d0)*sigma0)
xerfmax2=((xI2(ixP2)+&
half*dxI)-xCent(2))/(sqrt(2.d0)*sigma0)
factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-&
erfc(xerfmax2))/4.d0
EM(ixP1,ixP2)=EM(ixP1,ixP2)+fluxSubC*factor
enddo !ixP2
enddo !ixP1
endif !emit
enddo !iSubC3
endif !smallflux
enddo !iSubC2
enddo !iSubC1
enddo
enddo
enddo !ix
deallocate(flux)
end subroutine integrate_emission_spherical