subroutine integrate_spectra_cartesian(igrid,wL,dwLg,xS,dxSg,spectra,numWL,&
numXS,fl)
integer, intent(in) :: igrid,numWL,numXS
double precision, intent(in) :: wL(numWL),xS(numXS)
double precision, intent(in) :: dwLg,dxSg
double precision, intent(inout) :: spectra(numWL,numXS)
type(te_fluid), intent(in) :: fl
integer :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,ixImin1,&
ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,ix1,ix2,ix3,ixOnew,j
double precision, allocatable :: flux(:,:,:),v(:,:,:),pth(:,:,:),Te(:,:,&
:),rho(:,:,:)
double precision :: wlc,wlwd,res,dst_slit,xslit,arcsec
double precision :: vloc(1:3),xloc(1:3),dxloc(1:3),xIloc(1:2),&
dxIloc(1:2)
integer :: nSubC1,nSubC2,nSubC3,iSubC1,iSubC2,iSubC3,iwL,ixS,ixSmin,&
ixSmax,iwLmin,iwLmax,nwL
double precision :: slit_width,dxSubC1,dxSubC2,dxSubC3,xerfmin1,xerfmin2,&
xerfmin3,xerfmax1,xerfmax2,xerfmax3,fluxSubC
double precision :: xSubC(1:3),xCent(1:2)
integer :: mass
double precision :: logTe,lineCent
character (30) :: ion
double precision :: spaceRsl,wlRsl,sigma_PSF,wslit
double precision :: sigma_wl,sigma_xs,factor
if (SI_unit) then
arcsec=7.25d5/unit_length
else
arcsec=7.25d7/unit_length
endif
if (activate_unit_arcsec) then
xslit=location_slit*arcsec
else
xslit=location_slit
endif
call get_line_info(spectrum_wl,ion,mass,logTe,lineCent,spaceRsl,wlRsl,&
sigma_PSF,wslit)
ixOmin1=ixmlo1
ixOmin2=ixmlo2
ixOmin3=ixmlo3
ixOmax1=ixmhi1
ixOmax2=ixmhi2
ixOmax3=ixmhi3
ixImin1=ixglo1
ixImin2=ixglo2
ixImin3=ixglo3
ixImax1=ixghi1
ixImax2=ixghi2
ixImax3=ixghi3
allocate(flux(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
v(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
pth(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
Te(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
rho(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3))
! get local EUV flux and velocity
call get_EUV(spectrum_wl,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 fl%get_pthermal(ps(igrid)%w,ps(igrid)%x,ixImin1,ixImin2,ixImin3,&
ixImax1,ixImax2,ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
ixOmax3,pth)
call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixImin1,ixImin2,ixImin3,ixImax1,&
ixImax2,ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,rho)
call fl%get_var_Rfactor(ps(igrid)%w,ps(igrid)%x,ixImin1,ixImin2,ixImin3,&
ixImax1,ixImax2,ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
ixOmax3,Te)
Te(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)=pth(ixOmin1:ixOmax1,&
ixOmin2:ixOmax2,ixOmin3:ixOmax3)/(Te(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)*rho(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3))
do ix1=ixOmin1,ixOmax1
do ix2=ixOmin2,ixOmax2
do ix3=ixOmin3,ixOmax3
do j=1,3
vloc(j)=ps(igrid)%w(ix1,ix2,ix3,iw_mom(j))/rho(ix1,ix2,ix3)
enddo
call dot_product_loc(vloc,vec_LOS,res)
v(ix1,ix2,ix3)=res
enddo
enddo
enddo
deallocate(rho)
slit_width=wslit*arcsec
sigma_wl=sigma_PSF*dwLg
sigma_xs=sigma_PSF*dxSg
do ix1=ixOmin1,ixOmax1
do ix2=ixOmin2,ixOmax2
do ix3=ixOmin3,ixOmax3
if (flux(ix1,ix2,ix3)>smalldouble) then
xloc(1:3)=ps(igrid)%x(ix1,ix2,ix3,1:3)
dxloc(1:3)=ps(igrid)%dx(ix1,ix2,ix3,1:3)
call get_cor_image(xloc,xIloc)
call dot_product_loc(dxloc,vec_xI1,res)
dxIloc(1)=abs(res)
if (xIloc(1)>=xslit-half*(slit_width+dxIloc(1)) .and. &
xIloc(1)<=xslit+half*(slit_width+dxIloc(1))) then
nSubC1=1;nSubC2=1;nSubC3=1;
nSubC1=max(nSubC1,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
1)*abs(vec_xI1(1))/(slit_width/16.d0)))
nSubC2=max(nSubC2,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
2)*abs(vec_xI1(2))/(slit_width/16.d0)))
nSubC3=max(nSubC3,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
3)*abs(vec_xI1(3))/(slit_width/16.d0)));
nSubC1=max(nSubC1,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
1)*abs(vec_xI2(1))/(dxSg/4.d0)))
nSubC2=max(nSubC2,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
2)*abs(vec_xI2(2))/(dxSg/4.d0)))
nSubC3=max(nSubC3,ceiling(ps(igrid)%dx(ix1,ix2,ix3,&
3)*abs(vec_xI2(3))/(dxSg/4.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;
! local line center and line width
if (SI_unit) then
fluxSubC=flux(ix1,ix2,ix3&
)*dxSubC1*dxSubC2*dxSubC3*unit_length*1.d2/dxSg/dxSg !DN s-1,s-1,s-1
wlc=lineCent*(1.d0+v(ix1,ix2,ix3)*unit_velocity*1.d2/const_c)
else
fluxSubC=flux(ix1,ix2,ix3&
)*dxSubC1*dxSubC2*dxSubC3*unit_length/dxSg/dxSg !DN s-1,s-1,s-1
wlc=lineCent*(1.d0+v(ix1,ix2,ix3)*unit_velocity/const_c)
endif
wlwd=sqrt(kb_cgs*Te(ix1,ix2,ix3)*unit_temperature/(mass*mp_cgs))
wlwd=wlwd*lineCent/const_c
!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)=xloc(1)-half*dxloc(1)+(iSubC1-half)*dxSubC1
xSubC(2)=xloc(2)-half*dxloc(2)+(iSubC2-half)*dxSubC2
xSubC(3)=xloc(3)-half*dxloc(3)+(iSubC3-half)*dxSubC3;
call get_cor_image(xSubC,xCent)
dst_slit=abs(xCent(1)-xslit) ! space distance to slit center
if (dst_slit<=half*slit_width) then
ixS=floor((xCent(2)-(xS(1)-half*dxSg))/dxSg)+1
ixSmin=max(1,ixS-3)
ixSmax=min(ixS+3,numXS)
iwL=floor((wlc-(wL(1)-half*dwLg))/dwLg)+1
nwL=3*ceiling(wlwd/dwLg+1)
iwLmin=max(1,iwL-nwL)
iwLmax=min(iwL+nwL,numWL)
! calculate the contribution to nearby pixels
do iwL=iwLmin,iwLmax
do ixS=ixSmin,ixSmax
xerfmin1=(wL(iwL)-half*dwLg-wlc)/sqrt(2.d0*(sigma_wl**2+&
wlwd**2))
xerfmax1=(wL(iwL)+half*dwLg-wlc)/sqrt(2.d0*(sigma_wl**2+&
wlwd**2))
xerfmin2=(xS(ixS)-half*dxSg-&
xCent(2))/(sqrt(2.d0)*sigma_xs)
xerfmax2=(xS(ixS)+half*dxSg-&
xCent(2))/(sqrt(2.d0)*sigma_xs)
factor=(erfc(xerfmin1)-erfc(xerfmax1))*(erfc(xerfmin2)-&
erfc(xerfmax2))/4.d0
spectra(iwL,ixS)=spectra(iwL,ixS)+fluxSubC*factor
enddo
enddo
! nearby pixels
endif
enddo
enddo
enddo
endif
endif
enddo
enddo
enddo
deallocate(flux,v,pth,Te)
end subroutine integrate_spectra_cartesian