Merge pull request #230 from danielver02/kgk-kperp-norm-250716

Kgk kperp norm 250716

Author: danielver02

input.md

@@ -62,6 +62,14 @@ Selection for root finding method.
 **`numiter`**  
 Maximum number of iterations in secant method.
 
+**`kperp_norm`**   
+Choice of:  
+
+- True: Follow's Stix (10-57) normalization convention.  
+- False: Multiplies Stix (10-57) by $k_{\perp}^2 d_{ref}^2$.
+
+Depending on the user's choice of normalization, `D_threshold` needs to be adjusted to account for additional factors of $k_{\perp}^6 d_{ref}^6$.
+
 **`D_threshold`**  
 Minimum threshold for secant method.
 

src/ALPS_NHDS.f90

@@ -60,6 +60,7 @@ subroutine calc_chi(chi,chi_low,j,kz,kperp,x)
   !! Subroutine that calculates the susceptibility of species j based on NHDS.
   use alps_var, only : bMnmaxs, bMBessel_zeros, bMbetas, bMalphas,bMpdrifts
   use alps_var, only : ms, qs, ns
+  use alps_var, only : kperp_norm
   implicit none
 
   double complex, intent(out) :: chi(3,3)
@@ -193,7 +194,11 @@ subroutine calc_chi(chi,chi_low,j,kz,kperp,x)
   chi(2,3)=Y(2,3)/(ell*ell)
   chi(3,1)=Y(3,1)/(ell*ell)
   chi(3,2)=Y(3,2)/(ell*ell)
-  chi(3,3)=kperp*kperp*2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))+Y(3,3)/(ell*ell)
+  if (kperp_norm) then
+     chi(3,3)=2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))+Y(3,3)/(ell*ell)
+  else
+     chi(3,3)=kperp*kperp*2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))+Y(3,3)/(ell*ell)
+  endif
 
   n=0
   chi_low(1,1,n)=Y0(1,1)/(ell*ell)
@@ -204,7 +209,11 @@ subroutine calc_chi(chi,chi_low,j,kz,kperp,x)
   chi_low(2,3,n)=Y0(2,3)/(ell*ell)
   chi_low(3,1,n)=Y0(3,1)/(ell*ell)
   chi_low(3,2,n)=Y0(3,2)/(ell*ell)
-  chi_low(3,3,n)=Y0(3,3)/(ell*ell)+kperp*kperp*2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))
+  if (kperp_norm) then
+     chi_low(3,3,n)=Y0(3,3)/(ell*ell)+2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))
+  else
+     chi_low(3,3,n)=Y0(3,3)/(ell*ell)+kperp*kperp*2.d0*x*vdrift/(ell*ell*kz*vtherm*vtherm*bMalphas(j))
+  endif
 
   n=1
   chi_low(1,1,n)=Y1(1,1)/(ell*ell)
@@ -242,6 +251,7 @@ subroutine calc_ypsilon(Y,j,n,kz,kperp,x)
   !! Calculates the Y-tensor according to Stix for a bi-Maxwelling, using the NHDS calculation.
   use alps_var, only : bMbetas, bMalphas,bMpdrifts
   use alps_var, only : ms, qs, ns
+  use alps_var, only : kperp_norm
   implicit none
 
   double complex, parameter ::  uniti=(0.d0,1.d0)
@@ -338,28 +348,29 @@ subroutine calc_ypsilon(Y,j,n,kz,kperp,x)
   endif
 
 
-  ! The tensor in Stix's (10-57)
-  !Y(1,1)=1.d0*(n*n)*BInz*An/z
-  !Y(1,2)=-uniti*n*(BInz-dBInzdz)*An
-  !Y(1,3)=kperp*n*BInz*Bn/(Omega*z)
-  !Y(2,1)=uniti*n*(BInz-dBInzdz)*An
-  !Y(2,2)=(1.d0*(n*n)*BInz/z+2.d0*z*BInz-2.d0*z*dBInzdz)*An
-  !Y(2,3)=uniti*kperp*(BInz-dBInzdz)*Bn/Omega
-  !Y(3,1)=kperp*BInz*n*Bn/(Omega*z)
-  !Y(3,2)=-uniti*kperp*(BInz-dBInzdz)*Bn/Omega
-  !Y(3,3)=2.d0*(x-1.d0*n*Omega)*BInz*Bn/(kz*vtherm*vtherm*bMalphas(j))
-
-  ! The tensor in Stix's (10-57), multiplied by kperp^2 d_ref^2
-  Y(1,1)=1.d0*(n*n)*BInz*An/zp
-  Y(1,2)=-uniti*n*(BInz-dBInzdz)*An*kperp*kperp
-  Y(1,3)=kperp*n*BInz*Bn/(Omega*zp)
-  Y(2,1)=uniti*n*(BInz-dBInzdz)*An*kperp*kperp
-  Y(2,2)=(1.d0*(n*n)*BInz/zp+kperp*kperp*2.d0*z*BInz-kperp*kperp*2.d0*z*dBInzdz)*An
-  Y(2,3)=uniti*kperp*kperp*kperp*(BInz-dBInzdz)*Bn/Omega
-  Y(3,1)=kperp*BInz*n*Bn/(Omega*zp)
-  Y(3,2)=-uniti*kperp*kperp*kperp*(BInz-dBInzdz)*Bn/Omega
-  Y(3,3)=kperp*kperp*2.d0*(x-1.d0*n*Omega)*BInz*Bn/(kz*vtherm*vtherm*bMalphas(j))
-
+  if (kperp_norm) then
+     ! The tensor in Stix's (10-57)
+     Y(1,1)=1.d0*(n*n)*BInz*An/z
+     Y(1,2)=-uniti*n*(BInz-dBInzdz)*An
+     Y(1,3)=kperp*n*BInz*Bn/(Omega*z)
+     Y(2,1)=uniti*n*(BInz-dBInzdz)*An
+     Y(2,2)=(1.d0*(n*n)*BInz/z+2.d0*z*BInz-2.d0*z*dBInzdz)*An
+     Y(2,3)=uniti*kperp*(BInz-dBInzdz)*Bn/Omega
+     Y(3,1)=kperp*BInz*n*Bn/(Omega*z)
+     Y(3,2)=-uniti*kperp*(BInz-dBInzdz)*Bn/Omega
+     Y(3,3)=2.d0*(x-1.d0*n*Omega)*BInz*Bn/(kz*vtherm*vtherm*bMalphas(j))
+  else
+     ! The tensor in Stix's (10-57), multiplied by kperp^2 d_ref^2
+     Y(1,1)=1.d0*(n*n)*BInz*An/zp
+     Y(1,2)=-uniti*n*(BInz-dBInzdz)*An*kperp*kperp
+     Y(1,3)=kperp*n*BInz*Bn/(Omega*zp)
+     Y(2,1)=uniti*n*(BInz-dBInzdz)*An*kperp*kperp
+     Y(2,2)=(1.d0*(n*n)*BInz/zp+kperp*kperp*2.d0*z*BInz-kperp*kperp*2.d0*z*dBInzdz)*An
+     Y(2,3)=uniti*kperp*kperp*kperp*(BInz-dBInzdz)*Bn/Omega
+     Y(3,1)=kperp*BInz*n*Bn/(Omega*zp)
+     Y(3,2)=-uniti*kperp*kperp*kperp*(BInz-dBInzdz)*Bn/Omega
+     Y(3,3)=kperp*kperp*2.d0*(x-1.d0*n*Omega)*BInz*Bn/(kz*vtherm*vtherm*bMalphas(j))
+  endif
 
   end subroutine
 
@@ -369,6 +380,7 @@ subroutine calc_chi_cold(chi,j,kz,kperp,x)
   !! Subroutine that calculates the susceptibility of species j based on the cold-plasma dispersion relation based on the paper Verscharen & Chandran, ApJ 764, 88, 2013.
   use alps_var, only : bMbetas,bMpdrifts
   use alps_var, only : ms, qs, ns
+  use alps_var, only : kperp_norm
   implicit none
 
   double complex, intent(out) :: chi(3,3)
@@ -434,16 +446,21 @@ subroutine calc_chi_cold(chi,j,kz,kperp,x)
   dispM=uniti*(1.d0/(ell*ell))*kperp*vdrift*Omega/((x-kz*vdrift)**2-Omega**2)
 
 
-  chi(1,1)=(dispR+dispL)/2.d0
-  chi(1,2)=-uniti*(dispR-dispL)/2.d0
-  chi(1,3)=dispJ
-  chi(2,1)=uniti*(dispR-dispL)/2.d0
-  chi(2,2)=(dispR+dispL)/2.d0
-  chi(2,3)=dispM
-  chi(3,1)=dispJ
-  chi(3,2)=-dispM
-  chi(3,3)=dispP
-
+  if (kperp_norm) then
+     chi(1,1)=(dispR+dispL)/2.d0
+     chi(1,2)=-uniti*(dispR-dispL)/2.d0
+     chi(1,3)=dispJ
+     chi(2,1)=uniti*(dispR-dispL)/2.d0
+     chi(2,2)=(dispR+dispL)/2.d0
+     chi(2,3)=dispM
+     chi(3,1)=dispJ
+     chi(3,2)=-dispM
+     chi(3,3)=dispP
+  else
+     write(*,*) 'ERROR: check kperp norm for cold plasma distribution!!'
+     stop
+  endif
+     
   end subroutine
 
 

src/ALPS_com.f90

@@ -32,7 +32,7 @@ subroutine pass_instructions
     use alps_var, only : ns, qs, ms, wroots, kperp, kpar, bessel_zero
     use alps_var, only : wave, chi0, chi0_low, pp, df0, vA, pi
     use alps_var, only : secant_method, numiter, D_threshold, D_gap, D_prec
-    use alps_var, only : use_map
+    use alps_var, only : use_map, kperp_norm
     use alps_var, only : ni, nr, omi, omf, gami, gamf, loggridg, loggridw
     use alps_var, only : determine_minima, n_resonance_interval, positions_principal, Tlim
     use alps_var, only : n_scan, scan, scan_option, relativistic, logfit, usebM
@@ -65,6 +65,7 @@ subroutine pass_instructions
     call mpi_bcast(D_gap,1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierror)
     call mpi_bcast(Tlim,       1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierror)
     call mpi_bcast(use_map, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierror)
+    call mpi_bcast(kperp_norm, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierror)
     call mpi_bcast(n_scan, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierror)
     call mpi_bcast(scan_option, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierror)
     call mpi_bcast(secant_method,    1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierror)