subroutine get_image(qunit,datatype,fl)
! integrate emission flux along line of sight (LOS)
! in a 3D simulation box and get a 2D EUV image
use mod_global_parameters
use mod_constants
integer, intent(in) :: qunit
type(te_fluid), intent(in) :: fl
character(20), intent(in) :: datatype
integer :: ix1,ix2,ix3,numXI1,numXI2,numWI
double precision :: xImin1,xImax1,xImin2,xImax2,xIcent1,xIcent2,dxI
double precision, allocatable :: xI1(:),xI2(:),dxI1(:),dxI2(:)
double precision, allocatable :: wI(:,:,:),wIs(:,:,:),EM(:,:),WLB(:,:,:)
double precision :: vec_temp1(1:3),vec_temp2(1:3)
double precision :: vec_z(1:3),vec_cor(1:3),xI_cor(1:2)
double precision :: res,LOS_psi,r_max,r_loc
integer :: mass
character (30) :: ion
double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
double precision :: arcsec,RHESSI_rsl,LASCO_rsl,pixel,R_occult,smallflux
integer :: iigrid,igrid,i,j,numSI,iw
logical :: emit
if (coordinate==spherical) then
call init_vectors_spherical()
else
! cartesian
call init_vectors_cartesian()
endif
! calculate domain of the image
if (coordinate==spherical) then
xImin1=-abs(xprobmax1)
xImin2=-abs(xprobmax1)
xImax1=abs(xprobmax1)
xImax2=abs(xprobmax1)
else
! calculate domain of the image
do ix1=1,2
if (ix1==1) vec_cor(1)=xprobmin1
if (ix1==2) vec_cor(1)=xprobmax1
do ix2=1,2
if (ix2==1) vec_cor(2)=xprobmin2
if (ix2==2) vec_cor(2)=xprobmax2
do ix3=1,2
if (ix3==1) vec_cor(3)=xprobmin3
if (ix3==2) vec_cor(3)=xprobmax3
if (big_image) then
r_loc=(vec_cor(1)-x_origin(1))**2
r_loc=r_loc+(vec_cor(2)-x_origin(2))**2
r_loc=r_loc+(vec_cor(3)-x_origin(3))**2
r_loc=sqrt(r_loc)
if (ix1==1 .and. ix2==1 .and. ix3==1) then
r_max=r_loc
else
r_max=max(r_max,r_loc)
endif
else
call get_cor_image(vec_cor,xI_cor)
if (ix1==1 .and. ix2==1 .and. ix3==1) then
xImin1=xI_cor(1)
xImax1=xI_cor(1)
xImin2=xI_cor(2)
xImax2=xI_cor(2)
else
xImin1=min(xImin1,xI_cor(1))
xImax1=max(xImax1,xI_cor(1))
xImin2=min(xImin2,xI_cor(2))
xImax2=max(xImax2,xI_cor(2))
endif
endif
enddo
enddo
enddo
if (big_image) then
xImin1=-r_max
xImin2=-r_max
xImax1=r_max
xImax2=r_max
endif
endif
xIcent1=(xImin1+xImax1)/2.d0
xIcent2=(xImin2+xImax2)/2.d0
! tables for image
if (SI_unit) then
arcsec=7.25d5/unit_length
else
arcsec=7.25d7/unit_length
endif
if (datatype=='image_euv') then
call get_line_info(wavelength,ion,mass,logTe,lineCent,spaceRsl,wlRsl,&
sigma_PSF,wslit)
dxI=spaceRsl*arcsec ! intrument resolution of image
smallflux=smalldouble
else if (datatype=='image_sxr') then
RHESSI_rsl=2.3d0/instrument_resolution_factor
dxI=RHESSI_rsl*arcsec
smallflux=1.d-40
else if (datatype=='image_whitelight') then
if (whitelight_instrument=='LASCO/C1') then
LASCO_rsl=5.6d0/instrument_resolution_factor
R_occult=1.1d0
else if (whitelight_instrument=='LASCO/C2') then
LASCO_rsl=11.4d0/instrument_resolution_factor
R_occult=2.d0
else if (whitelight_instrument=='LASCO/C3') then
LASCO_rsl=56.d0/instrument_resolution_factor
R_occult=3.7d0
else
call MPISTOP('Whitelight synthesis: instrument is not supported!')
endif
dxI=LASCO_rsl*arcsec
if (R_occultor>1.d0) R_occult=R_occultor
R_occult=R_occult*const_Rsun/unit_length
smallflux=1.d-20
endif
numXI1=8*ceiling((xImax1-xIcent1)/dxI/8.d0)
xImin1=xIcent1-numXI1*dxI
xImax1=xIcent1+numXI1*dxI
numXI1=numXI1*2
numXI2=8*ceiling((xImax2-xIcent2)/dxI/8.d0)
xImin2=xIcent2-numXI2*dxI
xImax2=xIcent2+numXI2*dxI
numXI2=numXI2*2
allocate(xI1(numXI1),xI2(numXI2),dxI1(numXI1),dxI2(numXI2))
do ix1=1,numXI1
xI1(ix1)=xImin1+dxI*(ix1-half)
dxI1(ix1)=dxI
enddo
do ix2=1,numXI2
xI2(ix2)=xImin2+dxI*(ix2-half)
dxI2(ix2)=dxI
enddo
! calculate emission
if (datatype=='image_euv' .or. datatype=='image_sxr') then
numWI=1
allocate(wI(numXI1,numXI2,numWI),wIs(numXI1,numXI2,numWI),EM(numXI1,&
numXI2))
wI=zero
wIs=zero
EM=zero
if (coordinate==cartesian) then
do iigrid=1,igridstail; igrid=igrids(iigrid);
call integrate_emission_cartesian(igrid,numXI1,numXI2,xI1,xI2,dxI,&
fl,datatype,EM)
enddo
else
do iigrid=1,igridstail; igrid=igrids(iigrid);
call integrate_emission_spherical(igrid,numXI1,numXI2,xI1,xI2,dxI,&
fl,datatype,EM)
enddo
endif
do ix1=1,numXI1
do ix2=1,numXI2
if (EM(ix1,ix2)>smallflux) wIs(ix1,ix2,1)=EM(ix1,ix2)
enddo
enddo
numSI=numXI1*numXI2*numWI
call MPI_ALLREDUCE(wIs,wI,numSI,MPI_DOUBLE_PRECISION,MPI_SUM,icomm,&
ierrmpi)
if (activate_unit_arcsec) then
xI1=xI1/arcsec
dxI1=dxI1/arcsec
xI2=xI2/arcsec
dxI2=dxI2/arcsec
endif
call output_data(qunit,xI1,xI2,dxI1,dxI2,wI,numXI1,numXI2,numWI,&
datatype)
deallocate(wI,wIs,EM)
else if (datatype=='image_whitelight') then
numWI=2
allocate(wI(numXI1,numXI2,numWI),wIs(numXI1,numXI2,numWI),WLB(numXI1,&
numXI2,numWI))
wI=zero
wIs=zero
WLB=zero
if (coordinate==spherical) then
do iigrid=1,igridstail; igrid=igrids(iigrid);
call integrate_whitelight_spherical(igrid,numXI1,numXI2,numWI,xI1,&
xI2,dxI,fl,datatype,WLB)
enddo
endif
do ix1=1,numXI1
do ix2=1,numXI2
if (WLB(ix1,ix2,1)>smallflux) then
wIs(ix1,ix2,1)=WLB(ix1,ix2,1)
wIs(ix1,ix2,2)=WLB(ix1,ix2,2)
endif
enddo
enddo
numSI=numXI1*numXI2*numWI
call MPI_ALLREDUCE(wIs,wI,numSI,MPI_DOUBLE_PRECISION,MPI_SUM,icomm,&
ierrmpi)
if (activate_unit_arcsec) then
xI1=xI1/arcsec
dxI1=dxI1/arcsec
xI2=xI2/arcsec
dxI2=dxI2/arcsec
endif
call output_data(qunit,xI1,xI2,dxI1,dxI2,wI,numXI1,numXI2,numWI,&
datatype)
deallocate(wI,wIs,WLB)
endif
deallocate(xI1,xI2,dxI1,dxI2)
end subroutine get_image