Update to JETPLUME- Key changes: output for degeneracy calculation and updated formulas for loading output (small but critical error!)

.gitignore

@@ -19,3 +19,5 @@ docs_out
 __pycache__
 # Python virtual environments
 .venv
+.DS_Store
+.ipynb_checkpoints/

linfpclib/linfpcplot.py

@@ -94,8 +94,8 @@ def plotlinfpc_gyro(linfpcdata,filename='',zoomin=False,vlim=None,plotlog=False,
         print("Changing ylim of plot!!!")
         plt.gca().set_ylim(ylim[0],zlim[1])
     plt.title(plttitle)
-    plt.xlabel('$v_{||}/v_{ts}$')
-    plt.ylabel('$v_{\perp}/v_{ts}$')
+    plt.xlabel(r'$v_{||}/v_{ts}$')
+    plt.ylabel(r'$v_{\perp}/v_{ts}$')
     plt.colorbar()
     plt.grid()
 
@@ -220,9 +220,12 @@ def plotlinfpc_gyro_dist(linfpcdata,filename,plotkey,zoomin=False,vlim=None,plot
         print("Changing ylim of plot!!!")
         plt.gca().set_ylim(ylim[0],zlim[1])
     plt.title(plttitle)
-    plt.xlabel('$v_{||}/v_{ts}$')
-    plt.ylabel('$v_{\perp}/v_{ts}$')
+    plt.xlabel(r'$v_{||}/v_{ts}$')
+    plt.ylabel(r'$v_{\perp}/v_{ts}$')
     plt.colorbar()
+    _cmap = plt.get_cmap("PuOr")
+    zero_color = _cmap(0.5)   # mid-point of the colormap = value 0
+    plt.gca().set_facecolor(zero_color)
     plt.grid()
 
     if(setequal):
@@ -239,7 +242,7 @@ def plotlinfpc_gyro_dist(linfpcdata,filename,plotkey,zoomin=False,vlim=None,plot
 
     plt.show()
 
-def plot_9pan_cart(foldername,filenametag,flnm='',specnum='01',computeEner=False, scalevelocity=1):
+def plot_9pan_cart(foldername,filenametag,dataoverwritepar=None,dataoverwriteperp1=None,dataoverwriteperp2=None,flnm='',specnum='01',computeEner=False, scalevelocity=1):
     """
     Makes 3x3 plot of projections of FPC vel signature in cartesian coordinates
 
@@ -257,16 +260,20 @@ def plot_9pan_cart(foldername,filenametag,flnm='',specnum='01',computeEner=False
 
     from linfpclib.linfpc import loadlinfpccart
 
-    flnmread = foldername + filenametag +'.cparcart.specie'+specnum+'.mode01'
-    print('Reading: ',flnmread)
-    cartpar = loadlinfpccart(flnmread)
+    if(dataoverwritepar is None):
+        flnmread = foldername + filenametag +'.cparcart.specie'+specnum+'.mode01'
+        print('Reading: ',flnmread)
+        cartpar = loadlinfpccart(flnmread)
 
-    flnmread = foldername + filenametag +'.cperp1.specie'+specnum+'.mode01'
-    cartperp1 = loadlinfpccart(flnmread)
-
-    flnmread = foldername + filenametag +'.cperp2.specie'+specnum+'.mode01'
-    cartperp2 = loadlinfpccart(flnmread)
+        flnmread = foldername + filenametag +'.cperp1.specie'+specnum+'.mode01'
+        cartperp1 = loadlinfpccart(flnmread)
 
+        flnmread = foldername + filenametag +'.cperp2.specie'+specnum+'.mode01'
+        cartperp2 = loadlinfpccart(flnmread)
+    else:
+        cartpar = dataoverwritepar
+        cartperp1 = dataoverwriteperp1
+        cartperp2 = dataoverwriteperp2
 
     fig, axs = plt.subplots(3,3,figsize=(3*5,3*5),sharex=True)
 
@@ -326,7 +333,7 @@ def plot_9pan_cart(foldername,filenametag,flnm='',specnum='01',computeEner=False
 
     plt.show()
 
-def plot_fs1_re_im_cart(foldername,filenametag,flnm='',specnum='01', scalevelocity=1):
+def plot_fs1_re_im_cart(foldername,filenametag,dataoverwritedistfunccart=None,flnm='',specnum='01', scalevelocity=1):
     """
     Makes 2x3 plot of projections of fs1 in cartesian coordinates
 
@@ -344,10 +351,14 @@ def plot_fs1_re_im_cart(foldername,filenametag,flnm='',specnum='01', scaleveloci
 
     from linfpclib.linfpc import loadlinfpccart_dist
 
-    flnmreadimag = foldername + filenametag +'.dfs.imag.specie'+specnum+'.mode01'
-    flnmreadreal = foldername + filenametag +'.dfs.real.specie'+specnum+'.mode01'
+    if(dataoverwritedistfunccart is None):
+        flnmreadimag = foldername + filenametag +'.dfs.imag.specie'+specnum+'.mode01'
+        flnmreadreal = foldername + filenametag +'.dfs.real.specie'+specnum+'.mode01'
 
-    distfunccart = loadlinfpccart_dist(flnmreadreal,flnmreadimag)
+        distfunccart = loadlinfpccart_dist(flnmreadreal,flnmreadimag)
+
+    else:
+        distfunccart = dataoverwritedistfunccart
 
     fig, axs = plt.subplots(2,3,figsize=(3*5,2*5),sharex=True)
 
@@ -378,6 +389,10 @@ def plot_fs1_re_im_cart(foldername,filenametag,flnm='',specnum='01', scaleveloci
 
             axs[_j,_i].axis('equal')
 
+            _cmap = plt.get_cmap("PuOr")
+            zero_color = _cmap(0.5)   # mid-point of the colormap = value 0
+            axs[_j,_i].set_facecolor(zero_color)
+
             fig.colorbar(_tempim, ax=axs[_j,_i])
 
             _j = _j + 1
@@ -590,3 +605,15 @@ def sweep2dplot(sweep2d,xkey,ykey,zkey,xlabel,ylabel,zlabel,flnm = '', xlim=[],y
         plt.savefig(flnm,format='png',dpi=300)
 
     plt.show()
+
+def compute_2d_from_3d_hist(out):
+    #out is data from 3d load for hist/fs1 data
+
+    out2df1ronly = lfpc.reduce_3d_to_projections(out['fs1_r'], out['vx'], out['vy'], out['vz'], 're_f')
+    out2dfironly = lfpc.reduce_3d_to_projections(out['fs1_i'], out['vx'], out['vy'], out['vz'], 'im_f')
+    out2d = out2df1ronly.copy()
+    out2d['im_fvxvy'] = out2dfironly['im_fvxvy'].copy()
+    out2d['im_fvxvz'] = out2dfironly['im_fvxvz'].copy()
+    out2d['im_fvyvz'] = out2dfironly['im_fvyvz'].copy()
+
+    return out2d

output.md

@@ -110,19 +110,29 @@ where `j` ranges from 1 to `nspec`.
 The normalization follows Eqns. 33-37 in [Klein, K. G., Howes, G. G.,
 and Brown, C. R., 2025](https://iopscience.iop.org/article/10.3847/2515-5172/add1c2)
 
-If `heat` is set to true, the next set of columns will be the power absorption or emission from each component. If `low_n` is set to true, additional terms associated with Landau, Transit time, and Cyclotron heating will be output. If `eigen` is false, this data will start in the 7th column. If eigen is true, this data will start in the 18+8 `nspec`+1st column.
+If `heat` is set to true, the next set of columns will be the power absorption or emission from each component. If `low_n` is set to true, additional terms associated with Landau, Transit time, and Cyclotron heating will be output. If `eigen` is false, this data will start in the 7th column. If `eigen` is true, this data will start in the 18+8 `nspec`+1st column.
 
-If `low_n` is true, we have
+If `low_n` is true, the ordering is:
+  1) All total powers P_j for j=1..nspec, then
+  2) For each species j, the six sub-terms (yy, yz, zy, zz, n=0, n=±1) consecutively.
+
+Concretely (1-based column numbers; {!eigen}=1 if eigen is false, else 0):
+
+- P_j: 18 + (8 - 6{!eigen})*nspec + j
+
+Let afterP = 18 + (8 - 6{!eigen})*nspec + nspec.
+Then for species j in {1,..,nspec}:
+
+- P_j^yy:   afterP + 6*(j-1) + 1   (Transit Time Damping term 1)
+- P_j^yz:   afterP + 6*(j-1) + 2   (Transit Time Damping term 2)
+- P_j^zy:   afterP + 6*(j-1) + 3   (Landau Damping term 1)
+- P_j^zz:   afterP + 6*(j-1) + 4   (Landau Damping term 2)
+- P_j^{n=0}: afterP + 6*(j-1) + 5  (sum of Landau and Transit Time Damping)
+- P_j^{n=+1}: afterP + 6*(j-1) + 6 (n=±1 Cyclotron Damping)
+- P_j^{n=-1}: afterP + 7*(j-1) + 7 (n=-1 Cyclotron Damping)
+
+Here, {!eigen} (negation of the eigen boolean) equals 1 if eigen is false and 0 if eigen is true.
 
-- 18+(8-6{!eigen})`nspec`+j. $P_j$
-- 18+(9-6{!eigen})`nspec`+j. $P_j^{yy}$ (Transit Time Damping term 1).
-- 18+(10-6{!eigen})`nspec`+j. $P_j^{yz}$ (Transit Time Damping term 2).
-- 18+(11-6{!eigen})`nspec`+j. $P_j^{zy}$ (Landau Damping term 1).
-- 18+(12-6{!eigen})`nspec`+j. $P_j^{zz}$ (Landau Damping term 2).
-- 18+(13-6{!eigen})`nspec`+j. $P_j^{n=0}$ (sum of Landau and Transit Time Damping).
-- 18+(14-6{!eigen})`nspec`+j. $P_j^{n=+1}$ ($n=+ 1$ Cyclotron Damping).
-- 18+(15-6{!eigen})`nspec`+j. $P_j^{n=-1}$ ($n=- 1$ Cyclotron Damping).
-Here, {!eigen} (negation of eigen boolean) is equal to 1 if eigen is false and 0 if eigen is true.
 
 This will be followed by 6`nspec` columns containing the parameter lists $\mathcal{P}_j$ for each species or component. 
 - 19+`noutperspec` `nspec`+6(j-1)-12({!eigen}). $T_{ref,\parallel}/T_{j,\parallel}$.
@@ -131,7 +141,7 @@ This will be followed by 6`nspec` columns containing the parameter lists $\mathc
 - 22+`noutperspec` `nspec`+6(j-1)-12({!eigen}). $q_{ref}/q_{j}$.
 - 23+`noutperspec` `nspec`+6(j-1)-12({!eigen}). $n_{j}/n_{ref}$.
 - 24+`noutperspec` `nspec`+6(j-1)-12({!eigen}). $v_{j,drift}/v_{A,ref}$.
-Here, noutperspec = 0 is the number of additional outputs created by setting heating or eigen to true. If eigen and heating are false, then noutperspec = 7, if eigen is false and heating is true, then noutperspec = 8, and if eigen is true and if eigen and heating are true, then noutperspec=15. Note that {!eigen} (negation of eigen boolean) is equal to 1 if eigen is false and 0 if eigen is true.
+Here, noutperspec = 0 is the number of additional outputs created by setting heating or eigen to true. If eigen and heating are false, then noutperspec = 8, if eigen is false and heating is true, then noutperspec = 8, and if eigen is true and if eigen and heating are true, then noutperspec=16. Note that {!eigen} (negation of eigen boolean) is equal to 1 if eigen is false and 0 if eigen is true.
 
 For ease of reference, consider the below tables for column numbers for two and three component calculations.
 

src/disprels.f90

@@ -412,23 +412,79 @@ end subroutine map_search
   !-=-=-=-=-
   subroutine test_disp
     !! Testing routine for single evaluation of dispersion relation.   
-    use vars, only : wroots, nroot_max, writeOut
+    use vars, only : wroots, nroot_max, writeOut, nspec, scan
     use vars, only : kperp, kpar, vtp
+
     implicit none
+
+    !Eigenfunctions
+     complex, dimension(1:3)       :: ef
+     !! Electric Field eigenfluctuation.
+
+     complex, dimension(1:3)       :: bf
+     !! Magnetic Field eigenfluctuation.
+
+     complex, dimension(1:nspec)     :: ns
+     !! Density eigenfluctuation.
+
+     complex, dimension(1:3,1:nspec) :: Us
+     !! Velocity eigenfluctuation.
+
+     !Heating
+     real, dimension(1:nspec) :: Ps
+     !! Power into/out of species/component in one wave period.
+
+     real, dimension(1:4,1:nspec) :: Ps_split
+     !! Power into/out of species/components, broken into different contributions.
+     !! Deprecated.
+
+
+
 
     complex :: omega
     !! Complex Frequency input.
 
     complex :: D
     !! Output of dispersion relation function.
+
+    integer :: i,j
 
+
+
     omega= cmplx(wroots(1,1),wroots(2,1))
 
     D = disp(omega)
-    
+
     write(*,'(4es14.4)') omega,D*vtp**6.
     write(*,'(2es14.4)')kperp,kpar
 
+    call calc_eigen(omega,ef,bf,Us,ns,Ps,Ps_split,.true.,.true.)
+
+
+    ! --- Print diagnostics ---
+    write(*,*) 'ef ='
+    write(*,'(3(1x,es14.6))') ef
+
+    write(*,*) 'bf ='
+    write(*,'(3(1x,es14.6))') bf
+
+    write(*,*) 'ns ='
+    write(*,'(100(1x,es14.6))') (ns(i), i=1,nspec)
+
+    write(*,*) 'Us ='
+    do i=1,nspec
+       write(*,'(3(1x,es14.6))') (Us(j,i), j=1,3)
+    end do
+
+    write(*,*) 'Ps ='
+    write(*,'(100(1x,es14.6))') (Ps(i), i=1,nspec)
+
+    write(*,*) 'Ps_split ='
+    do i=1,nspec
+       write(*,'(4(1x,es14.6))') (Ps_split(j,i), j=1,4)
+    end do
+
+
   end subroutine test_disp
 
 !-=-=-=-=-
@@ -3378,31 +3434,36 @@ subroutine get_double_out_name(outName,tensorName,fmt,fmt_tnsr,out_type,diff)
      if (scan(2)%eigen_s) then
         if (scan(2)%heat_s) then
            if (low_n) then
-              !>>>GGH: 1/18/23
-                 write(fmt,'(a,i0,a)')'(6es15.6e3,12es15.6e3,',21*nspec,'es15.6e3)'
+              ! eigen + heat + low_n  (MATCH single: 22*nspec)
+              write(fmt,'(a,i0,a)') '(6es15.6e3,12es15.6e3,', 22*nspec, 'es15.6e3)'
            else
-              write(fmt,'(a,i0,a)')'(6es15.6e3,12es15.6e3,',15*nspec,'es15.6e3,es15.6e3)'
+              ! eigen + heat (non-low_n)
+              write(fmt,'(a,i0,a)') '(6es15.6e3,12es15.6e3,', 15*nspec, 'es15.6e3,es15.6e3)'
            endif
-           out_type=0
+           out_type = 0
         else
-           write(fmt,'(a,i0,a)')'(6es15.6e3,12es15.6e3,',14*nspec,'es15.6e3)'
-           out_type=1
+           ! eigen only
+           write(fmt,'(a,i0,a)') '(6es15.6e3,12es15.6e3,', 14*nspec, 'es15.6e3)'
+           out_type = 1
         endif
      else
         if (scan(2)%heat_s) then
            if (low_n) then
-              !>>>GGH: 1/18/23
-                 write(fmt,'(a,i0,a)')'(6es15.6e3,',13*nspec,'es15.6e3)'
+              ! heat only + low_n  (MATCH single: 14*nspec)
+              write(fmt,'(a,i0,a)') '(6es15.6e3,', 14*nspec, 'es15.6e3)'
            else
-              write(fmt,'(a,i0,a)')'(6es15.6e3,',7*nspec,'es15.6e3,es15.6e3)'
+              ! heat only (non-low_n)  (remove extra trailing es field)
+              write(fmt,'(a,i0,a)') '(6es15.6e3,', 7*nspec, 'es15.6e3)'
            endif
-           out_type=2
+           out_type = 2
         else
-           write(fmt,'(a,i0,a)')'(6es15.6e3,',6*nspec,'es15.6e3)'
-           out_type=3
+          ! frequency only
+          write(fmt,'(a,i0,a)') '(6es15.6e3,', 6*nspec, 'es15.6e3)'
+          out_type = 3
         endif
      endif
 
+
      if ((scan(1)%tensor_s).and.(scan(2)%tensor_s)) &
           write(fmt_tnsr,'(a,i0,a)')'(4es15.6e3,',18*nspec,'es15.6e3)'