subroutine get_EUV(wl,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,w,x,fl,flux)
! calculate the local emission intensity of given EUV line (optically thin)
! wavelength is the wave length of the emission line
! unit [DN cm^-1 s^-1 pixel^-1]
! ingrate flux along line of sight: DN s^-1 pixel^-1
use mod_global_parameters
integer, intent(in) :: wl
integer, intent(in) :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
double precision, intent(in) :: x(ixImin1:ixImax1,ixImin2:ixImax2,&
ixImin3:ixImax3,1:ndim)
double precision, intent(in) :: w(ixImin1:ixImax1,ixImin2:ixImax2,&
ixImin3:ixImax3,1:nw)
type(te_fluid), intent(in) :: fl
double precision, intent(out) :: flux(ixImin1:ixImax1,ixImin2:ixImax2,&
ixImin3:ixImax3)
integer :: n_table
double precision, allocatable :: t_table(:),f_table(:)
integer :: ix1,ix2,ix3,iTt,i
double precision :: pth(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
Te(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3),&
Ne(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3)
double precision :: logT,logGT,GT
! selecting emission table
select case(wl)
case(94)
n_table=n_aia
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_aia(1:n_aia)
f_table(1:n_table)=f_94(1:n_aia)
case(131)
n_table=n_aia
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_aia(1:n_aia)
f_table(1:n_table)=f_131(1:n_aia)
case(171)
n_table=n_aia
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_aia(1:n_aia)
f_table(1:n_table)=f_171(1:n_aia)
case(193)
n_table=n_aia
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_aia(1:n_aia)
f_table(1:n_table)=f_193(1:n_aia)
case(211)
n_table=n_aia
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_aia(1:n_aia)
f_table(1:n_table)=f_211(1:n_aia)
case(304)
n_table=n_aia
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_aia(1:n_aia)
f_table(1:n_table)=f_304(1:n_aia)
case(335)
n_table=n_aia
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_aia(1:n_aia)
f_table(1:n_table)=f_335(1:n_aia)
case(1354)
n_table=n_iris
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_iris(1:n_iris)
f_table(1:n_table)=f_1354(1:n_iris)
case(263)
n_table=n_eis
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_eis1(1:n_eis)
f_table(1:n_table)=f_263(1:n_eis)
case(264)
n_table=n_eis
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_eis2(1:n_eis)
f_table(1:n_table)=f_264(1:n_eis)
case(192)
n_table=n_eis
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_eis2(1:n_eis)
f_table(1:n_table)=f_192(1:n_eis)
case(255)
n_table=n_eis
allocate(t_table(1:n_table))
allocate(f_table(1:n_table))
t_table(1:n_table)=t_eis2(1:n_eis)
f_table(1:n_table)=f_255(1:n_eis)
case default
allocate(t_table(1))
allocate(f_table(1))
call mpistop("Unknown wavelength")
end select
call fl%get_pthermal(w,x,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,pth)
call fl%get_rho(w,x,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,Ne)
call fl%get_var_Rfactor(w,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)/(Ne(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)*Te(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3))*unit_temperature
if (SI_unit) then
Ne(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)=Ne(ixOmin1:ixOmax1,&
ixOmin2:ixOmax2,ixOmin3:ixOmax3)*unit_numberdensity/1.d6 !m-3,m-3,m-3 -> cm-3
flux(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)=Ne(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)**2
else
Ne(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3)=Ne(ixOmin1:ixOmax1,&
ixOmin2:ixOmax2,ixOmin3:ixOmax3)*unit_numberdensity
flux(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)=Ne(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3)**2
endif
select case(wl)
case(94,131,171,193,211,304,335,1354)
! temperature table in log
do i=1,n_table
if(f_table(i) .gt. 1.d-99) then
f_table(i)=log10(f_table(i))
else
f_table(i)=-99.d0
endif
enddo
logGT=zero
do ix3=ixOmin3,ixOmax3
do ix2=ixOmin2,ixOmax2
do ix1=ixOmin1,ixOmax1
logT=log10(Te(ix1,ix2,ix3))
if (logT>=t_table(1) .and. logT<=t_table(n_table)) then
do iTt=1,n_table-1
if (logT>=t_table(iTt) .and. logT<t_table(iTt+1)) then
logGT=f_table(iTt)*(logT-t_table(iTt+&
1))/(t_table(iTt)-t_table(iTt+1))+f_table(iTt+&
1)*(logT-t_table(iTt))/(t_table(iTt+1)-t_table(iTt))
endif
enddo
flux(ix1,ix2,ix3)=flux(ix1,ix2,ix3)*(10**(logGT))
if(flux(ix1,ix2,ix3)<smalldouble) flux(ix1,ix2,ix3)=0.d0
else
flux(ix1,ix2,ix3)=zero
endif
enddo
enddo
enddo
case default
! temperature table linear
do ix3=ixOmin3,ixOmax3
do ix2=ixOmin2,ixOmax2
do ix1=ixOmin1,ixOmax1
if (Te(ix1,ix2,ix3)>=t_table(1) .and. Te(ix1,ix2,&
ix3)<=t_table(n_table)) then
do iTt=1,n_table-1
if (Te(ix1,ix2,ix3)>=t_table(iTt) .and. Te(ix1,ix2,&
ix3)<t_table(iTt+1)) then
GT=f_table(iTt)*(Te(ix1,ix2,&
ix3)-t_table(iTt+1))/(t_table(iTt)-t_table(iTt+1))+&
f_table(iTt+1)*(Te(ix1,ix2,&
ix3)-t_table(iTt))/(t_table(iTt+1)-t_table(iTt))
endif
enddo
flux(ix1,ix2,ix3)=flux(ix1,ix2,ix3)*GT
if(flux(ix1,ix2,ix3)<smalldouble) flux(ix1,ix2,ix3)=0.d0
else
flux(ix1,ix2,ix3)=zero
endif
enddo
enddo
enddo
end select
deallocate(t_table,f_table)
end subroutine get_EUV