SUBROUTINE angular_density(params,spp)
IMPLICIT NONE
TYPE(KORC_PARAMS), INTENT(IN) :: params
TYPE(SPECIES), DIMENSION(:), ALLOCATABLE, INTENT(IN) :: spp
CHARACTER(MAX_STRING_LENGTH) :: var_name
LOGICAL, DIMENSION(:,:,:), ALLOCATABLE :: bool_pixel_array
LOGICAL :: bool
REAL(rp), DIMENSION(3) :: binorm, n, nperp
REAL(rp), DIMENSION(3) :: X, V, B, E, XC
REAL(rp), DIMENSION(:,:,:), ALLOCATABLE :: angle_pixel_array
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: np_angular_pixel
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: Psyn_angular_pixel
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: np_lambda_pixel
REAL(rp), DIMENSION(:,:,:,:), ALLOCATABLE :: Psyn_lambda_pixel
REAL(rp) :: q, m, k, u, g, l, threshold_angle
REAL(rp) :: psi, chi, beta, theta
REAL(rp), DIMENSION(:), ALLOCATABLE :: P_lambda
REAL(rp) :: P_angular
REAL(rp), DIMENSION(:), ALLOCATABLE :: photon_energy
REAL(rp) :: r,solid_angle
REAL(rp) :: angle, clockwise
REAL(rp) :: units
REAL(rp), DIMENSION(:), ALLOCATABLE :: Psyn_send_buffer,Psyn_receive_buffer, np_send_buffer, np_receive_buffer
REAL(rp) :: lc
REAL(rp) :: dSA ! Element of solid angle
REAL(rp), DIMENSION(:), ALLOCATABLE :: zeta
INTEGER :: ii,jj,ll,ss,pp
INTEGER :: numel, mpierr
ALLOCATE(bool_pixel_array(cam%num_pixels(1),cam%num_pixels(2),2))
! (NX,NY,2)
ALLOCATE(angle_pixel_array(cam%num_pixels(1),cam%num_pixels(2),2))
! (NX,NY,2)
ALLOCATE(np_angular_pixel(cam%num_pixels(1),cam%num_pixels(2), &
cam%Nlambda,params%num_species))
ALLOCATE(Psyn_angular_pixel(cam%num_pixels(1),cam%num_pixels(2), &
cam%Nlambda,params%num_species))
ALLOCATE(np_lambda_pixel(cam%num_pixels(1),cam%num_pixels(2), &
cam%Nlambda,params%num_species))
ALLOCATE(Psyn_lambda_pixel(cam%num_pixels(1),cam%num_pixels(2), &
cam%Nlambda,params%num_species))
ALLOCATE(P_lambda(cam%Nlambda))
ALLOCATE(zeta(cam%Nlambda))
ALLOCATE(photon_energy(cam%Nlambda))
np_angular_pixel = 0.0_rp
Psyn_angular_pixel = 0.0_rp
np_lambda_pixel = 0.0_rp
Psyn_lambda_pixel = 0.0_rp
P_lambda = 0.0_rp
zeta = 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,n,nperp,X,XC,V,B,E, &
!$OMP& bool_pixel_array,angle_pixel_array,k,u,g,l,threshold_angle,theta,&
!$OMP& psi,chi,beta,P_lambda,bool,angle,clockwise,ii,jj, &
!$OMP& ll,pp,r,lc,zeta,solid_angle,dSA)&
!$OMP& SHARED(params,spp,ss,Psyn_angular_pixel,np_angular_pixel, &
!$OMP& np_lambda_pixel,Psyn_lambda_pixel)
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) ! Critical wavelength
zeta = lc/cam%lambda
binorm = binorm/SQRT(DOT_PRODUCT(binorm,binorm))
threshold_angle = (1.5_rp*k*cam%lambda_max/C_PI)**(1.0_rp/3.0_rp)
! In radians
call check_if_visible(X,V/u,threshold_angle,bool,angle)
if (bool.EQV..TRUE.) then
X(1:2) = clockwise_rotation(X(1:2),angle)
V(1:2) = clockwise_rotation(V(1:2),angle)
binorm(1:2) = clockwise_rotation(binorm(1:2),angle)
call calculate_rotation_angles(X,bool_pixel_array, &
angle_pixel_array)
clockwise = ATAN2(X(2),X(1))
if (clockwise.LT.0.0_rp) clockwise = clockwise + 2.0_rp*C_PI
do ii=1_idef,cam%num_pixels(1) ! NX
do jj=1_idef,cam%num_pixels(2) ! NY
if (bool_pixel_array(ii,jj,1)) then
angle = angle_pixel_array(ii,jj,1) - clockwise
XC = (/cam%position(1)*COS(angle),-cam%position(1)* &
SIN(angle),cam%position(2)/)
n = XC - X
r = SQRT(DOT_PRODUCT(n,n))
n = n/r
dSA = DOT_PRODUCT(n,XC/SQRT(DOT_PRODUCT(XC,XC)))* &
cam%pixel_area/r**2
beta = ACOS(DOT_PRODUCT(n,binorm))
if (beta.GT.0.5_rp*C_PI) psi = beta - 0.5_rp*C_PI
if (beta.LT.0.5_rp*C_PI) psi = 0.5_rp*C_PI - beta
nperp = n - DOT_PRODUCT(n,binorm)*binorm
nperp = nperp/SQRT(DOT_PRODUCT(nperp,nperp))
chi = ABS(ACOS(DOT_PRODUCT(nperp,V/u)))
theta = ABS(ACOS(DOT_PRODUCT(n,V/u)))
if (theta .LE. threshold_angle) then
call P_integral(zeta,P_lambda)
P_lambda = (C_C*C_E**2)*P_lambda/(SQRT(3.0_rp)* &
C_E0*g**2*cam%lambda**3)
P_lambda = dSA*P_lambda/(2.0_rp*C_PI*(1.0_rp - &
COS(1.0_rp/g)))
if (cam%photon_count) then
Psyn_lambda_pixel(ii,jj,:,ss) = &
Psyn_lambda_pixel(ii,jj,:,ss) + &
P_lambda/photon_energy
else
Psyn_lambda_pixel(ii,jj,:,ss) = &
Psyn_lambda_pixel(ii,jj,:,ss) + P_lambda
end if
np_lambda_pixel(ii,jj,:,ss) = &
np_lambda_pixel(ii,jj,:,ss) + 1.0_rp
end if
do ll=1_idef,cam%Nlambda ! Nlambda
if ((chi.LT.chic(g,k,cam%lambda(ll))).AND. &
(psi.LT.psic(k,cam%lambda(ll)))) then
P_angular = Psyn(g,psi,k,cam%lambda(ll),chi)
if (P_angular.GT.0.0_rp) then
P_angular = dSA*P_angular
if (cam%photon_count) then
Psyn_angular_pixel(ii,jj,ll,ss) = &
Psyn_angular_pixel(ii,jj,ll,ss) &
+ P_angular/photon_energy(ll)
else
Psyn_angular_pixel(ii,jj,ll,ss) = &
Psyn_angular_pixel(ii,jj,ll,ss) + &
P_angular
end if
np_angular_pixel(ii,jj,ll,ss) = &
np_angular_pixel(ii,jj,ll,ss) + 1.0_rp
end if
end if
end do ! Nlambda
end if
if (bool_pixel_array(ii,jj,2)) then
angle = angle_pixel_array(ii,jj,2) - clockwise
XC = (/cam%position(1)*COS(angle),-cam%position(1)* &
SIN(angle),cam%position(2)/)
n = XC - X
r = SQRT(DOT_PRODUCT(n,n))
n = n/r
dSA = DOT_PRODUCT(n,XC/SQRT(DOT_PRODUCT(XC,XC)))* &
cam%pixel_area/r**2
beta = ACOS(DOT_PRODUCT(n,binorm))
if (beta.GT.0.5_rp*C_PI) psi = beta - 0.5_rp*C_PI
if (beta.LT.0.5_rp*C_PI) psi = 0.5_rp*C_PI - beta
nperp = n - DOT_PRODUCT(n,binorm)*binorm
nperp = nperp/SQRT(DOT_PRODUCT(nperp,nperp))
chi = ABS(ACOS(DOT_PRODUCT(nperp,V/u)))
theta = ABS(ACOS(DOT_PRODUCT(n,V/u)))
if (theta .LE. threshold_angle) then
call P_integral(zeta,P_lambda)
P_lambda = (C_C*C_E**2)*P_lambda/(SQRT(3.0_rp)* &
C_E0*g**2*cam%lambda**3)
P_lambda = dSA*P_lambda/(2.0_rp*C_PI*(1.0_rp - &
COS(1.0_rp/g)))
if (cam%photon_count) then
Psyn_lambda_pixel(ii,jj,:,ss) = &
Psyn_lambda_pixel(ii,jj,:,ss) + &
P_lambda/photon_energy
else
Psyn_lambda_pixel(ii,jj,:,ss) = &
Psyn_lambda_pixel(ii,jj,:,ss) + P_lambda
end if
np_lambda_pixel(ii,jj,:,ss) = &
np_lambda_pixel(ii,jj,:,ss) + 1.0_rp
end if
do ll=1_idef,cam%Nlambda ! Nlambda
if ((chi.LT.chic(g,k,cam%lambda(ll))).AND. &
(psi.LT.psic(k,cam%lambda(ll)))) then
P_angular = Psyn(g,psi,k,cam%lambda(ll),chi)
if (P_angular.GT.0.0_rp) then
P_angular = dSA*P_angular
if (cam%photon_count) then
Psyn_angular_pixel(ii,jj,ll,ss) = &
Psyn_angular_pixel(ii,jj,ll,ss) &
+ P_angular/photon_energy(ll)
else
Psyn_angular_pixel(ii,jj,ll,ss) = &
Psyn_angular_pixel(ii,jj,ll,ss) + &
P_angular
end if
np_angular_pixel(ii,jj,ll,ss) = &
np_angular_pixel(ii,jj,ll,ss) + 1.0_rp
end if
end if
end do ! Nlambda
end if
end do ! NY
end do ! NX
end if ! check if bool == TRUE
end if ! if confined
end do ! particles
!$OMP END PARALLEL DO
end do ! species
! * * * * * * * IMPORTANT * * * * * * *
! * * * * * * * * * * * * * * * * * * *
! Here Psyn has units of Watts/m
! or (photons/s)(m^-1). See above.
! * * * * * * * * * * * * * * * * * * *
! * * * * * * * IMPORTANT * * * * * * *
units = 1.0_rp
if (params%mpi_params%nmpi.GT.1_idef) then
numel = cam%num_pixels(1)*cam%num_pixels(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_angular_pixel,(/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_angular_pixel,(/numel/))
CALL MPI_REDUCE(np_send_buffer,np_receive_buffer,numel, &
MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
Psyn_angular_pixel = RESHAPE(Psyn_receive_buffer, &
(/cam%num_pixels(1),cam%num_pixels(2),cam%Nlambda, &
params%num_species/))
np_angular_pixel = RESHAPE(np_receive_buffer, &
(/cam%num_pixels(1),cam%num_pixels(2),cam%Nlambda, &
params%num_species/))
var_name = 'np_angular_pixel'
call save_snapshot_var(params,np_angular_pixel,var_name)
var_name = 'Psyn_angular_pixel'
call save_snapshot_var(params,Psyn_angular_pixel,var_name)
end if
Psyn_send_buffer = RESHAPE(Psyn_lambda_pixel,(/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_lambda_pixel,(/numel/))
CALL MPI_REDUCE(np_send_buffer,np_receive_buffer,numel, &
MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
Psyn_lambda_pixel = RESHAPE(Psyn_receive_buffer, &
(/cam%num_pixels(1),cam%num_pixels(2),cam%Nlambda, &
params%num_species/))
np_lambda_pixel = RESHAPE(np_receive_buffer, &
(/cam%num_pixels(1),cam%num_pixels(2),cam%Nlambda, &
params%num_species/))
var_name = 'np_lambda_pixel'
call save_snapshot_var(params,np_lambda_pixel,var_name)
var_name = 'Psyn_lambda_pixel'
call save_snapshot_var(params,Psyn_lambda_pixel,var_name)
end if
DEALLOCATE(Psyn_send_buffer)
DEALLOCATE(Psyn_receive_buffer)
DEALLOCATE(np_send_buffer)
DEALLOCATE(np_receive_buffer)
CALL MPI_BARRIER(MPI_COMM_WORLD,mpierr)
else
var_name = 'np_angular_pixel'
call save_snapshot_var(params,np_angular_pixel,var_name)
var_name = 'np_lambda_pixel'
call save_snapshot_var(params,np_lambda_pixel,var_name)
var_name = 'Psyn_angular_pixel'
Psyn_angular_pixel = units*Psyn_angular_pixel
call save_snapshot_var(params,Psyn_angular_pixel,var_name)
var_name = 'Psyn_lambda_pixel'
call save_snapshot_var(params,Psyn_lambda_pixel,var_name)
end if
DEALLOCATE(bool_pixel_array)
DEALLOCATE(angle_pixel_array)
DEALLOCATE(np_angular_pixel)
DEALLOCATE(Psyn_angular_pixel)
DEALLOCATE(np_lambda_pixel)
DEALLOCATE(Psyn_lambda_pixel)
DEALLOCATE(P_lambda)
DEALLOCATE(zeta)
DEALLOCATE(photon_energy)
END SUBROUTINE angular_density