SUBROUTINE integrated_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), DIMENSION(:,:), ALLOCATABLE :: P_lambda_pixel
REAL(rp), DIMENSION(:,:), ALLOCATABLE :: P_angular_pixel
REAL(rp), DIMENSION(:,:), ALLOCATABLE :: P_l_pixel
REAL(rp), DIMENSION(:,:), ALLOCATABLE :: P_a_pixel
REAL(rp), DIMENSION(:), ALLOCATABLE :: np_pixel
REAL(rp), DIMENSION(:), ALLOCATABLE :: P_lambda, P_angular
REAL(rp) :: q, m, k, u, g, l, threshold_angle, threshold_angle_simple_model
REAL(rp) :: psi, chi, beta, theta
REAL(rp), DIMENSION(:), ALLOCATABLE :: photon_energy
REAL(rp) :: r,solid_angle
REAL(rp) :: angle, clockwise
REAL(rp) :: units
REAL(rp), DIMENSION(:), ALLOCATABLE :: send_buffer, 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), &
params%num_species))
ALLOCATE(Psyn_angular_pixel(cam%num_pixels(1),cam%num_pixels(2), &
params%num_species))
ALLOCATE(np_lambda_pixel(cam%num_pixels(1),cam%num_pixels(2), &
params%num_species))
ALLOCATE(Psyn_lambda_pixel(cam%num_pixels(1),cam%num_pixels(2), &
params%num_species))
ALLOCATE(P_lambda(cam%Nlambda))
ALLOCATE(P_angular(cam%Nlambda))
ALLOCATE(P_l_pixel(cam%Nlambda,params%num_species))
ALLOCATE(P_a_pixel(cam%Nlambda,params%num_species))
ALLOCATE(np_pixel(params%num_species))
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_l_pixel = 0.0_rp
P_a_pixel = 0.0_rp
np_pixel = 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, &
!$OMP& threshold_angle_simple_model,theta,&
!$OMP& psi,chi,beta,bool,angle,clockwise,ii,jj,ll,pp,r,lc,zeta, &
!$OMP& P_lambda,P_angular,solid_angle,dSA)&
!$OMP& SHARED(params,spp,ss,Psyn_angular_pixel,np_angular_pixel, &
!$OMP& np_lambda_pixel,Psyn_lambda_pixel,P_l_pixel,P_a_pixel,np_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
threshold_angle_simple_model = 1.0_rp/g
np_pixel(ss) = np_pixel(ss) + 1.0_rp ! We count all the confined particles.
call check_if_visible(X,V/u,MAXVAL((/threshold_angle,threshold_angle_simple_model/)),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 ! Here, angle is modified w.r.t. check_if_visible.
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)))
P_lambda = 0.0_rp
P_angular = 0.0_rp
if (theta .LE. threshold_angle_simple_model) 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)
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(ll) = Psyn(g,psi,k,cam%lambda(ll),chi)
if (P_angular(ll).GT.0.0_rp) then
np_angular_pixel(ii,jj,ss) = np_angular_pixel(ii,jj,ss) + 1.0_rp
else
P_angular(ll) = 0.0_rp
end if
end if
end do ! Nlambda
P_lambda = dSA*P_lambda/(2.0_rp*C_PI*(1.0_rp - COS(1.0_rp/g)))
P_angular = dSA*P_angular
P_l_pixel(:,ss) = P_l_pixel(:,ss) + P_lambda
P_a_pixel(:,ss) = P_a_pixel(:,ss) + P_angular
if (cam%photon_count) then
P_lambda = P_lambda/photon_energy
P_angular = P_angular/photon_energy
Psyn_lambda_pixel(ii,jj,ss) = Psyn_lambda_pixel(ii,jj,ss) + SUM(P_lambda)
Psyn_angular_pixel(ii,jj,ss) = Psyn_angular_pixel(ii,jj,ss) + SUM(P_angular)
else
Psyn_lambda_pixel(ii,jj,ss) = Psyn_lambda_pixel(ii,jj,ss) + trapz(cam%lambda,P_lambda)
Psyn_angular_pixel(ii,jj,ss) = Psyn_angular_pixel(ii,jj,ss) + trapz(cam%lambda,P_angular)
end if
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)))
P_lambda = 0.0_rp
P_angular = 0.0_rp
if (theta .LE. threshold_angle_simple_model) 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)
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(ll) = Psyn(g,psi,k,cam%lambda(ll),chi)
if (P_angular(ll).GT.0.0_rp) then
np_angular_pixel(ii,jj,ss) = np_angular_pixel(ii,jj,ss) + 1.0_rp
else
P_angular(ll) = 0.0_rp
end if
end if
end do ! Nlambda
P_lambda = dSA*P_lambda/(2.0_rp*C_PI*(1.0_rp - COS(1.0_rp/g)))
P_angular = dSA*P_angular
P_l_pixel(:,ss) = P_l_pixel(:,ss) + P_lambda
P_a_pixel(:,ss) = P_a_pixel(:,ss) + P_angular
if (cam%photon_count) then
P_lambda = P_lambda/photon_energy
P_angular = P_angular/photon_energy
Psyn_lambda_pixel(ii,jj,ss) = Psyn_lambda_pixel(ii,jj,ss) + SUM(P_lambda)
Psyn_angular_pixel(ii,jj,ss) = Psyn_angular_pixel(ii,jj,ss) + SUM(P_angular)
else
Psyn_lambda_pixel(ii,jj,ss) = Psyn_lambda_pixel(ii,jj,ss) + trapz(cam%lambda,P_lambda)
Psyn_angular_pixel(ii,jj,ss) = Psyn_angular_pixel(ii,jj,ss) + trapz(cam%lambda,P_angular)
end if
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)*params%num_species
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(Psyn_angular_pixel,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,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(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),params%num_species/))
var_name = 'Psyn_angular_pixel'
call save_snapshot_var(params,Psyn_angular_pixel,var_name)
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(np_angular_pixel,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
np_angular_pixel = RESHAPE(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),params%num_species/))
var_name = 'np_angular_pixel'
call save_snapshot_var(params,np_angular_pixel,var_name)
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(Psyn_lambda_pixel,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,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(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),params%num_species/))
var_name = 'Psyn_lambda_pixel'
call save_snapshot_var(params,Psyn_lambda_pixel,var_name)
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(np_lambda_pixel,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
np_lambda_pixel = RESHAPE(receive_buffer,(/cam%num_pixels(1),cam%num_pixels(2),params%num_species/))
var_name = 'np_lambda_pixel'
call save_snapshot_var(params,np_lambda_pixel,var_name)
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
numel = params%num_species
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = np_pixel
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
np_pixel = receive_buffer
var_name = 'np_pixel'
call save_snapshot_var(params,np_pixel,var_name)
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
numel = cam%Nlambda*params%num_species
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(P_a_pixel,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
P_a_pixel = RESHAPE(receive_buffer,(/cam%Nlambda,params%num_species/))
var_name = 'P_a_pixel'
call save_snapshot_var(params,P_a_pixel,var_name)
end if
DEALLOCATE(send_buffer)
DEALLOCATE(receive_buffer)
ALLOCATE(send_buffer(numel))
ALLOCATE(receive_buffer(numel))
send_buffer = RESHAPE(P_l_pixel,(/numel/))
receive_buffer = 0.0_rp
CALL MPI_REDUCE(send_buffer,receive_buffer,numel,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
P_l_pixel = RESHAPE(receive_buffer,(/cam%Nlambda,params%num_species/))
var_name = 'P_l_pixel'
call save_snapshot_var(params,P_l_pixel,var_name)
end if
DEALLOCATE(send_buffer)
DEALLOCATE(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'
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)
var_name = 'np_pixel'
call save_snapshot_var(params,np_pixel,var_name)
var_name = 'P_a_pixel'
call save_snapshot_var(params,P_a_pixel,var_name)
var_name = 'P_l_pixel'
call save_snapshot_var(params,P_l_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(P_angular)
DEALLOCATE(P_l_pixel)
DEALLOCATE(P_a_pixel)
DEALLOCATE(np_pixel)
DEALLOCATE(zeta)
DEALLOCATE(photon_energy)
END SUBROUTINE integrated_angular_density