subroutine integrate_emission_cartesian(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 :: xbmin1,xbmin2,xbmin3,xbmax1,xbmax2,xbmax3,xd1,xd2,&
xd3
double precision, allocatable :: flux(:,:,:)
double precision :: res
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
double precision :: xSubC(1:3),dxSubC1,dxSubC2,dxSubC3,xCent(1:2)
integer :: mass
double precision :: logTe
character (30) :: ion
double precision :: lineCent
double precision :: sigma_PSF,spaceRsl,wlRsl,sigma0,factor,wslit
double precision :: arcsec,pixel,RHESSI_rsl,area_1AU
double precision :: aa,bb
ixOmin1=ixmlo1
ixOmin2=ixmlo2
ixOmin3=ixmlo3
ixOmax1=ixmhi1
ixOmax2=ixmhi2
ixOmax3=ixmhi3
ixImin1=ixglo1
ixImin2=ixglo2
ixImin3=ixglo3
ixImax1=ixghi1
ixImax2=ixghi2
ixImax3=ixghi3
xbmin1=rnode(rpxmin1_,igrid)
xbmin2=rnode(rpxmin2_,igrid)
xbmin3=rnode(rpxmin3_,igrid)
xbmax1=rnode(rpxmax1_,igrid)
xbmax2=rnode(rpxmax2_,igrid)
xbmax3=rnode(rpxmax3_,igrid)
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
do ix1=ixOmin1,ixOmax1
do ix2=ixOmin2,ixOmax2
do ix3=ixOmin3,ixOmax3
nSubC1=1;nSubC2=1;nSubC3=1;
nSubC1=max(nSubC1,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
1)*abs(vec_xI1(1))/(dxI/2.d0)))
nSubC2=max(nSubC2,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
2)*abs(vec_xI1(2))/(dxI/2.d0)))
nSubC3=max(nSubC3,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
3)*abs(vec_xI1(3))/(dxI/2.d0)));
nSubC1=max(nSubC1,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
1)*abs(vec_xI2(1))/(dxI/2.d0)))
nSubC2=max(nSubC2,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
2)*abs(vec_xI2(2))/(dxI/2.d0)))
nSubC3=max(nSubC3,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
3)*abs(vec_xI2(3))/(dxI/2.d0)));
dxSubC1=ps(igrid)%dx(ix1,ix2,ix3,1)/nSubC1
dxSubC2=ps(igrid)%dx(ix1,ix2,ix3,2)/nSubC2
dxSubC3=ps(igrid)%dx(ix1,ix2,ix3,3)/nSubC3;
if (datatype=='image_euv') then
if (SI_unit) then
fluxSubC=flux(ix1,ix2,ix3&
)*dxSubC1*dxSubC2*dxSubC3*unit_length*1.d2/dxI/dxI !DN s-1,s-1,s-1
else
fluxSubC=flux(ix1,ix2,ix3&
)*dxSubC1*dxSubC2*dxSubC3*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)*dxSubC1*dxSubC2*dxSubC3*unit_length**&
3/area_1AU
endif
if (fluxSubC>smalldouble) then
!dividing a cell to several parts to get more accurate integrating values
do iSubC1=1,nSubC1
do iSubC2=1,nSubC2
do iSubC3=1,nSubC3
xSubC(1)=ps(igrid)%x(ix1,ix2,ix3,1)-half*ps(igrid)%dx(ix1,ix2,ix3,&
1)+(iSubC1-half)*dxSubC1
xSubC(2)=ps(igrid)%x(ix1,ix2,ix3,2)-half*ps(igrid)%dx(ix1,ix2,ix3,&
2)+(iSubC2-half)*dxSubC2
xSubC(3)=ps(igrid)%x(ix1,ix2,ix3,3)-half*ps(igrid)%dx(ix1,ix2,ix3,&
3)+(iSubC3-half)*dxSubC3;
! mapping the 3D coordinate to location at the image
call get_cor_image(xSubC,xCent)
! 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
enddo
enddo
enddo !iSubC
endif
enddo
enddo
enddo !ix
deallocate(flux)
end subroutine integrate_emission_cartesian