subroutine integrate_whitelight_spherical(igrid,numXI1,numXI2,numWI,xI1,&
xI2,dxI,fl,datatype,WLB)
integer, intent(in) :: igrid,numXI1,numXI2,numWI
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) :: WLB(numXI1,numXI2,numWI)
integer :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,ixO1,ixO2,&
ixO3,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,ix1,ix2,ix3,i,j
double precision, allocatable :: flux(:,:,:),Ne(:,:,:)
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,RsubC
double precision :: sigma_PSF,sigma0,arcsec,pixel,LASCO_rsl
double precision :: A,B,C,D,Rc,Ne0,TBsubC,PBsubC,factor
double precision :: R_thick,dotp,dvolume,R_occult
double precision :: xSubC(1:3),dxSubC1,dxSubC2,dxSubC3,xCent(1:2),&
xSubC_car(1:3)
double precision :: dxl(1:3),x_sph(1:3),dx_sph(1:3)
double precision :: unitv_r(1:3),unitv_theta(1:3),unitv_phi(1:3)
logical :: emit
ixOmin1=ixmlo1;ixOmin2=ixmlo2;ixOmin3=ixmlo3;
ixOmax1=ixmhi1;ixOmax2=ixmhi2;ixOmax3=ixmhi3;
ixImin1=ixglo1;ixImin2=ixglo2;ixImin3=ixglo3;
ixImax1=ixghi1;ixImax2=ixghi2;ixImax3=ixghi3;
if (SI_unit) then
arcsec=7.25d5/unit_length
else
arcsec=7.25d7/unit_length
endif
allocate(Ne(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3))
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
endif
if (R_occultor>1.d0) R_occult=R_occultor
R_occult=R_occult*const_Rsun/unit_length
call fl%get_rho(ps(igrid)%w,ps(igrid)%x,ixImin1,ixImin2,ixImin3,ixImax1,&
ixImax2,ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,Ne)
sigma_PSF=1.d0
pixel=LASCO_rsl*arcsec
sigma0=sigma_PSF*pixel
! integrate emission
R_thick=R_opt_thick*const_Rsun/unit_length
do ix1=ixOmin1,ixOmax1
do ix2=ixOmin2,ixOmax2
do ix3=ixOmin3,ixOmax3
x_sph(1:3)=ps(igrid)%x(ix1,ix2,ix3,1:3)
dx_sph(1:3)=ps(igrid)%dx(ix1,ix2,ix3,1:3)
dxl(1)=dx_sph(1) ! cell size in length
dxl(2)=x_sph(1)*dx_sph(2) ! cell size in length
dxl(3)=x_sph(1)*dsin(x_sph(2))*dx_sph(3) ! cell size in length
Ne0=Ne(ix1,ix2,ix3)*unit_numberdensity
!dividing a cell to several sub-cells to get more accurate integrating values
nSubC1=1;nSubC2=1;nSubC3=1;
call get_unit_vector_spherical(x_sph,unitv_r,unitv_theta,unitv_phi)
call dot_product_loc(unitv_r,vec_xI1,dotp)
nSubC1=max(nSubC1,ceiling(dxl(1)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_r,vec_xI2,dotp)
nSubC1=max(nSubC1,ceiling(dxl(1)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_theta,vec_xI1,dotp)
nSubC2=max(nSubC2,ceiling(dxl(2)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_theta,vec_xI2,dotp)
nSubC2=max(nSubC2,ceiling(dxl(2)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_phi,vec_xI1,dotp)
nSubC3=max(nSubC3,ceiling(dxl(3)*abs(dotp)/(dxI/2.d0)))
call dot_product_loc(unitv_phi,vec_xI2,dotp)
nSubC3=max(nSubC3,ceiling(dxl(3)*abs(dotp)/(dxI/2.d0)))
! integrate sub-cells
do iSubc1=1,nSubC1
! sub-cell center coordinate in spherical
xSubC(1)=x_sph(1)-half*dx_sph(1)+(iSubC1-half)*dx_sph(1)/nSubC1
RsubC=xSubC(1)
dxSubC1=dx_sph(1)/nSubC1 ! sub-cell size in length
call get_Thomson_parameters(RsubC,A,B,C,D)
do iSubc2=1,nSubC2
! sub-cell center coordinate in spherical
xSubC(2)=x_sph(2)-half*dx_sph(2)+(iSubC2-half)*dx_sph(2)/nSubC2
dxSubC2=xSubC(1)*dx_sph(2)/nSubC2 ! sub-cell size in length
dxSubC3=xSubC(1)*dsin(xSubC(2))*dx_sph(3)/nSubC3 !sub-cell size in length
dvolume=dxSubC1*dxSubC2*dxSubC3
do iSubc3=1,nSubC3
! sub-cell center coordinate in spherical
xSubC(3)=x_sph(3)-half*dx_sph(3)+(iSubC3-half)*dx_sph(3)/nSubC3
call get_cor_image_spherical(xSubC,xCent)
Rc=dsqrt(xCent(1)**2+xCent(2)**2) !distance to sun center (on the image plane)
! whether the local emitted photons can arrive the telescope
emit=.false.
if (Rc>R_occult) then
emit=.true.
! scaterring flux from cm-3,cm-3,cm-3 of plasma
call get_whitelight_Thomson(RsubC,Rc,Ne0,A,B,C,D,TBsubC,&
PBsubC)
TBsubC=TBsubC*dvolume*unit_length/dxI/dxI
PBsubC=PBsubC*dvolume*unit_length/dxI/dxI
if (TBsubC<1.d-20) emit=.false.
endif
if (emit) then
! mapping the 3D coordinate to location at the image
! 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
WLB(ixP1,ixP2,1)=WLB(ixP1,ixP2,1)+TBsubC*factor
WLB(ixP1,ixP2,2)=WLB(ixP1,ixP2,2)+PBsubC*factor
enddo !ixP2
enddo !ixP1
endif
enddo !iSubC3
enddo !iSubC2
enddo !iSubC1
enddo
enddo
enddo !ix
deallocate(Ne)
end subroutine integrate_whitelight_spherical