subroutine get_spectrum(qunit,datatype,fl)
integer, intent(in) :: qunit
character(20), intent(in) :: datatype
type(te_fluid), intent(in) :: fl
integer :: numWL,numXS,iwL,ixS,numWI,ix1,ix2,ix3
double precision :: dwLg,dxSg,xSmin,xSmax,xScent,wLmin,wLmax
double precision, allocatable :: wL(:),xS(:),dwL(:),dxS(:)
double precision, allocatable :: wI(:,:,:),spectra(:,:),spectra_rc(:,:)
double precision :: vec_cor(1:3),xI_cor(1:2)
double precision :: res,r_loc,r_max
integer :: mass
character (30) :: ion
double precision :: logTe,lineCent,sigma_PSF,spaceRsl,wlRsl,wslit
double precision :: unitv,arcsec,RHESSI_rsl,pixel
integer :: iigrid,igrid,i,j,numS
double precision :: xLmin,xLmax,xslit
if (coordinate==spherical) then
call init_vectors_spherical()
else
! cartesian
call init_vectors_cartesian()
endif
! calculate domain in space
if (coordinate==spherical) then
xSmin=-abs(xprobmax1)
xSmax=abs(xprobmax1)
else
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
xSmin=xI_cor(2)
xSmax=xI_cor(2)
else
xSmin=min(xSmin,xI_cor(2))
xSmax=max(xSmax,xI_cor(2))
endif
endif
enddo
enddo
enddo
if (big_image) then
xSmin=-r_max
xSmax=r_max
endif
endif
xScent=(xSmin+xSmax)/2.d0
! tables for storing spectra data
if (SI_unit) then
arcsec=7.25d5/unit_length
else
arcsec=7.25d7/unit_length
endif
call get_line_info(spectrum_wl,ion,mass,logTe,lineCent,spaceRsl,wlRsl,&
sigma_PSF,wslit)
dxSg=spaceRsl*arcsec
numXS=ceiling((xSmax-xScent)/dxSg)
xSmin=xScent-numXS*dxSg
xSmax=xScent+numXS*dxSg
numXS=numXS*2
dwLg=wlRsl
numWL=2*int((spectrum_window_max-spectrum_window_min)/(2.d0*dwLg))
wLmin=(spectrum_window_max+spectrum_window_min)/2.d0-dwLg*numWL/2
wLmax=(spectrum_window_max+spectrum_window_min)/2.d0+dwLg*numWL/2
allocate(wL(numWL),dwL(numWL),xS(numXS),dxS(numXS))
numWI=1
allocate(wI(numWL,numXS,numWI),spectra(numWL,numXS),spectra_rc(numWL,&
numXS))
do iWL=1,numWL
wL(iwL)=wLmin+iwL*dwLg-half*dwLg
dwL=dwLg
enddo
do ixS=1,numXS
xS(ixS)=xSmin+dxSg*(ixS-half)
dxS(ixS)=dxSg
enddo
! find slit and do integration
spectra=zero
do iigrid=1,igridstail; igrid=igrids(iigrid);
do ix1=1,2
if (ix1==1) vec_cor(1)=rnode(rpxmin1_,igrid)
if (ix1==2) vec_cor(1)=rnode(rpxmax1_,igrid)
do ix2=1,2
if (ix2==1) vec_cor(2)=rnode(rpxmin2_,igrid)
if (ix2==2) vec_cor(2)=rnode(rpxmax2_,igrid)
do ix3=1,2
if (ix3==1) vec_cor(3)=rnode(rpxmin3_,igrid)
if (ix3==2) vec_cor(3)=rnode(rpxmax3_,igrid)
call get_cor_image(vec_cor,xI_cor)
if (ix1==1 .and. ix2==1 .and. ix3==1) then
xLmin=xI_cor(1)
xLmax=xI_cor(1)
else
xLmin=min(xLmin,xI_cor(1))
xLmax=max(xLmax,xI_cor(1))
endif
enddo
enddo
enddo
if (activate_unit_arcsec) then
xslit=location_slit*arcsec
else
xslit=location_slit
endif
if (xslit>=xLmin-wslit*arcsec .and. xslit<=xLmax+wslit*arcsec) then
call integrate_spectra_cartesian(igrid,wL,dwLg,xS,dxSg,spectra,numWL,&
numXS,fl)
endif
enddo
numS=numWL*numXS
call MPI_ALLREDUCE(spectra,spectra_rc,numS,MPI_DOUBLE_PRECISION, MPI_SUM,&
icomm,ierrmpi)
do iwL=1,numWL
do ixS=1,numXS
if (spectra_rc(iwL,ixS)>smalldouble) then
wI(iwL,ixS,1)=spectra_rc(iwL,ixS)
else
wI(iwL,ixS,1)=zero
endif
enddo
enddo
if (activate_unit_arcsec) then
xS=xS/arcsec
dxS=dxS/arcsec
endif
call output_data(qunit,wL,xS,dwL,dxS,wI,numWL,numXS,numWI,datatype)
deallocate(wL,xS,dwL,dxS,spectra,spectra_rc,wI)
end subroutine get_spectrum