SUBROUTINE spectral_density(params,spp)
IMPLICIT NONE
TYPE(KORC_PARAMS), INTENT(IN) :: params
TYPE(SPECIES), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: spp
CHARACTER(MAX_STRING_LENGTH) :: var_name
REAL(rp), DIMENSION(3) :: binorm
REAL(rp), DIMENSION(3) :: X, V, B, E
REAL(rp), DIMENSION(:), ALLOCATABLE :: P
REAL(rp), DIMENSION(:), ALLOCATABLE :: zeta
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: np
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: Psyn_lambda
REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: PTot
REAL(rp) :: Rpol, Zpol
REAL(rp) :: q, m, k, u, g, lc
REAL(rp), DIMENSION(:), ALLOCATABLE :: photon_energy
INTEGER :: ii,jj,ll,ss,pp
REAL(rp), DIMENSION(:), ALLOCATABLE :: Psyn_send_buffer,Psyn_receive_buffer
REAL(rp), DIMENSION(:), ALLOCATABLE :: np_send_buffer,np_receive_buffer
REAL(rp), DIMENSION(:), ALLOCATABLE :: PTot_send_buffer,PTot_receive_buffer
INTEGER :: numel, mpierr
REAL(rp) :: units
ALLOCATE(np(pplane%grid_dims(1),pplane%grid_dims(2),cam%Nlambda,params%num_species))
ALLOCATE(Psyn_lambda(pplane%grid_dims(1),pplane%grid_dims(2),cam%Nlambda,params%num_species))
ALLOCATE(PTot(pplane%grid_dims(1),pplane%grid_dims(2),params%num_species))
ALLOCATE(P(cam%Nlambda))
ALLOCATE(zeta(cam%Nlambda))
ALLOCATE(photon_energy(cam%Nlambda))
np = 0.0_rp
Psyn_lambda = 0.0_rp
P = 0.0_rp
PTot = 0.0_rp
photon_energy = C_h*C_C/cam%lambda
do ss=1_idef,params%num_species
q = ABS(spp(ss)%q)*params%cpp%charge
m = spp(ss)%m*params%cpp%mass
!$OMP PARALLEL DO FIRSTPRIVATE(q,m,photon_energy) &
!$OMP& PRIVATE(binorm,X,V,B,E,k,u,g,lc,ii,jj,ll,pp,zeta,P,Rpol,Zpol) &
!$OMP& SHARED(params,spp,ss,Psyn_lambda,PTot,np)
do pp=1_idef,spp(ss)%ppp
if ( spp(ss)%vars%flag(pp) .EQ. 1_idef ) then
V = spp(ss)%vars%V(pp,:)*params%cpp%velocity
X = spp(ss)%vars%X(pp,:)*params%cpp%length
g = spp(ss)%vars%g(pp)
B = spp(ss)%vars%B(pp,:)*params%cpp%Bo
E = spp(ss)%vars%E(pp,:)*params%cpp%Eo
binorm = cross(V,E) + cross(V,cross(V,B))
u = SQRT(DOT_PRODUCT(V,V))
k = q*SQRT(DOT_PRODUCT(binorm,binorm))/(spp(ss)%vars%g(pp)*m*u**3)
lc = (4.0_rp*C_PI/3.0_rp)/(k*g**3)
zeta = lc/cam%lambda
call P_integral(zeta,P)
P = (C_C*C_E**2)*P/(SQRT(3.0_rp)*C_E0*g**2*cam%lambda**3)
Rpol = SQRT(SUM(X(1:2)**2))
Zpol = X(3)
ii = FLOOR((Rpol - pplane%Rmin)/pplane%DR) + 1_idef
jj = FLOOR((Zpol + ABS(pplane%Zmin))/pplane%DZ) + 1_idef
Psyn_lambda(ii,jj,:,ss) = Psyn_lambda(ii,jj,:,ss) + P
np(ii,jj,:,ss) = np(ii,jj,:,ss) + 1_idef
PTot(ii,jj,ss) = PTot(ii,jj,ss) + spp(ss)%vars%Prad(pp);
end if ! if confined
end do ! particles
!$OMP END PARALLEL DO
end do ! species
! * * * * * * * IMPORTANT * * * * * * *
! * * * * * * * * * * * * * * * * * * *
! Here Psyn_lambda has units of Watts/m
! or (photons/s)(m^-1). See above.
! * * * * * * * * * * * * * * * * * * *
! * * * * * * * IMPORTANT * * * * * * *
if (params%mpi_params%nmpi.GT.1_idef) then
numel = pplane%grid_dims(1)*pplane%grid_dims(2)*cam%Nlambda*params%num_species
ALLOCATE(Psyn_send_buffer(numel))
ALLOCATE(Psyn_receive_buffer(numel))
ALLOCATE(np_send_buffer(numel))
ALLOCATE(np_receive_buffer(numel))
Psyn_send_buffer = RESHAPE(Psyn_lambda,(/numel/))
CALL MPI_REDUCE(Psyn_send_buffer,Psyn_receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
np_send_buffer = RESHAPE(np,(/numel/))
CALL MPI_REDUCE(np_send_buffer,np_receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
numel = pplane%grid_dims(1)*pplane%grid_dims(2)*params%num_species
ALLOCATE(PTot_send_buffer(numel))
ALLOCATE(PTot_receive_buffer(numel))
PTot_send_buffer = RESHAPE(PTot,(/numel/))
CALL MPI_REDUCE(PTot_send_buffer,PTot_receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
Psyn_lambda = RESHAPE(Psyn_receive_buffer,(/pplane%grid_dims(1),pplane%grid_dims(2),cam%Nlambda,params%num_species/))
np = RESHAPE(np_receive_buffer,(/pplane%grid_dims(1),pplane%grid_dims(2),cam%Nlambda,params%num_species/))
var_name = 'np_pplane'
call save_snapshot_var(params,np,var_name)
var_name = 'Psyn_pplane'
call save_snapshot_var(params,Psyn_lambda,var_name)
PTot = RESHAPE(PTot_receive_buffer,(/pplane%grid_dims(1),pplane%grid_dims(2),params%num_species/))
units = params%cpp%mass*(params%cpp%velocity**3)/params%cpp%length
PTot = units*PTot ! (Watts)
var_name = 'PTot_pplane'
call save_snapshot_var(params,PTot,var_name)
end if
DEALLOCATE(Psyn_send_buffer)
DEALLOCATE(Psyn_receive_buffer)
DEALLOCATE(PTot_send_buffer)
DEALLOCATE(PTot_receive_buffer)
DEALLOCATE(np_send_buffer)
DEALLOCATE(np_receive_buffer)
CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
else
var_name = 'np_pplane'
call save_snapshot_var(params,np,var_name)
var_name = 'Psyn_pplane'
call save_snapshot_var(params,Psyn_lambda,var_name)
units = params%cpp%mass*(params%cpp%velocity**3)/params%cpp%length
PTot = units*PTot ! (Watts)
var_name = 'PTot_pplane'
call save_snapshot_var(params,PTot,var_name)
end if
DEALLOCATE(np)
DEALLOCATE(Psyn_lambda)
DEALLOCATE(PTot)
DEALLOCATE(P)
DEALLOCATE(zeta)
DEALLOCATE(photon_energy)
END SUBROUTINE spectral_density