From 87e09d5133be4c48266fe832d75db26fbfb253f0 Mon Sep 17 00:00:00 2001 From: kgklein Date: Wed, 17 Sep 2025 16:37:32 -0700 Subject: [PATCH] Updated heating output to differentiate between n=1 and n=-1, rather than combining them in one output column. This changes the column structure. See output.md for details. --- data/PLUME-Output-Table-2.svg | 8 +- data/PLUME-Output-Table.svg | 8 +- output.md | 5 +- plotter/example/example_kpar.plt | 61 ++++++--- plotter/example/example_kperp.plt | 59 ++++++--- src/disprels.f90 | 207 +++++++++++++++++++----------- src/fpc.f90 | 45 ++++--- src/functions.f90 | 2 +- src/vars.f90 | 2 +- 9 files changed, 248 insertions(+), 149 deletions(-) diff --git a/data/PLUME-Output-Table-2.svg b/data/PLUME-Output-Table-2.svg index c9e4b0a..b98f457 100644 --- a/data/PLUME-Output-Table-2.svg +++ b/data/PLUME-Output-Table-2.svg @@ -1,14 +1,14 @@ - + - + - + - + diff --git a/data/PLUME-Output-Table.svg b/data/PLUME-Output-Table.svg index 9dd985e..d47579b 100644 --- a/data/PLUME-Output-Table.svg +++ b/data/PLUME-Output-Table.svg @@ -1,14 +1,14 @@ - + - + - + - + diff --git a/output.md b/output.md index 09ee23a..0594e49 100644 --- a/output.md +++ b/output.md @@ -120,7 +120,8 @@ If `low_n` is true, we have - 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=\pm 1}$ ($n=\pm 1$ Cyclotron 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. @@ -132,6 +133,8 @@ This will be followed by 6`nspec` columns containing the parameter lists $\mathc - 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. +For ease of reference, consider the below tables for column numbers for two and three component calculations. + ![Table for output column numbering for a two-component plasma.](./data/PLUME-Output-Table.svg) ![Table for output column numbering for a three-component plasma.](./data/PLUME-Output-Table-2.svg) diff --git a/plotter/example/example_kpar.plt b/plotter/example/example_kpar.plt index 466a7d6..36584c4 100644 --- a/plotter/example/example_kpar.plt +++ b/plotter/example/example_kpar.plt @@ -43,6 +43,7 @@ dy=0.225 base_color='#000000' n0_color='#009e73' n1_color='#e68f00' +nm1_color='#00f86e' set log x set xrange [k_min:k_max] @@ -135,7 +136,9 @@ file_in(mm) u 2:($35/sgn($5)) w l lc rgb base_color dt 2 title '$\gamma_p<0$' ,\ file_in(mm) u 2:(-$41/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ file_in(mm) u 2:($41/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ file_in(mm) u 2:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 2:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 2:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 2:(-$43/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 2:($43/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Electron Damping Rates unset label @@ -151,10 +154,12 @@ set label '$|\gamma_e|/|\omega_{\textrm{r}}|$' at graph -0.3,0.15 rotate by 90 plot \ file_in(mm) u 2:(-$36/sgn($5)) w l lc rgb base_color dt 1 title '$\gamma_e>0$' ,\ file_in(mm) u 2:($36/sgn($5)) w l lc rgb base_color dt 2 title '$\gamma_e<0$' ,\ -file_in(mm) u 2:(-$47/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ -file_in(mm) u 2:($47/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ -file_in(mm) u 2:(-$48/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 2:($48/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 2:(-$48/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ +file_in(mm) u 2:($48/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ +file_in(mm) u 2:(-$49/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ +file_in(mm) u 2:($49/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 2:(-$50/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 2:($50/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Polarization unset label @@ -224,7 +229,9 @@ file_in(mm) u 2:($35/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ file_in(mm) u 2:(-$41/sgn($5)) w l lc rgb n0_color dt 1 title '$\gamma_p^{n=0}>0$' ,\ file_in(mm) u 2:($41/sgn($5)) w l lc rgb n0_color dt 2 title '$\gamma_p^{n=0}<0$' ,\ file_in(mm) u 2:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 2:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 2:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 2:(-$43/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 2:($43/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Electron Damping Rates unset label @@ -239,10 +246,12 @@ set yrange [gam_min:gam_max] plot \ file_in(mm) u 2:(-$36/sgn($5)) w l lc rgb base_color dt 1 title '' ,\ file_in(mm) u 2:($36/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ -file_in(mm) u 2:(-$47/sgn($5)) w l lc rgb n0_color dt 1 title '$\gamma_e^{n=0}>0$' ,\ -file_in(mm) u 2:($47/sgn($5)) w l lc rgb n0_color dt 2 title '$\gamma_e^{n=0}<0$' ,\ -file_in(mm) u 2:(-$48/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 2:($48/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 2:(-$48/sgn($5)) w l lc rgb n0_color dt 1 title '$\gamma_e^{n=0}>0$' ,\ +file_in(mm) u 2:($48/sgn($5)) w l lc rgb n0_color dt 2 title '$\gamma_e^{n=0}<0$' ,\ +file_in(mm) u 2:(-$49/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ +file_in(mm) u 2:($49/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 2:(-$50/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 2:($50/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Polarization unset label @@ -310,8 +319,10 @@ file_in(mm) u 2:(-$35/sgn($5)) w l lc rgb base_color dt 1 title '' ,\ file_in(mm) u 2:($35/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ file_in(mm) u 2:(-$41/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ file_in(mm) u 2:($41/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ -file_in(mm) u 2:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '$\gamma_p^{n=\pm 1}>0$' ,\ -file_in(mm) u 2:($42/sgn($5)) w l lc rgb n1_color dt 2 title '$\gamma_p^{n=\pm 1}<0$' +file_in(mm) u 2:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '$\gamma_p^{n=+ 1}>0$' ,\ +file_in(mm) u 2:($42/sgn($5)) w l lc rgb n1_color dt 2 title '$\gamma_p^{n=+ 1}<0$' ,\ +file_in(mm) u 2:(-$43/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 2:($43/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Electron Damping Rates unset label @@ -326,10 +337,12 @@ set yrange [gam_min:gam_max] plot \ file_in(mm) u 2:(-$36/sgn($5)) w l lc rgb base_color dt 1 title '' ,\ file_in(mm) u 2:($36/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ -file_in(mm) u 2:(-$47/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ -file_in(mm) u 2:($47/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ -file_in(mm) u 2:(-$48/sgn($5)) w l lc rgb n1_color dt 1 title '$\gamma_e^{n=\pm 1}>0$' ,\ -file_in(mm) u 2:($48/sgn($5)) w l lc rgb n1_color dt 2 title '$\gamma_e^{n=\pm 1}<0$' +file_in(mm) u 2:(-$48/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ +file_in(mm) u 2:($48/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ +file_in(mm) u 2:(-$49/sgn($5)) w l lc rgb n1_color dt 1 title '$\gamma_e^{n=+ 1}>0$' ,\ +file_in(mm) u 2:($49/sgn($5)) w l lc rgb n1_color dt 2 title '$\gamma_e^{n=+ 1}<0$' ,\ +file_in(mm) u 2:(-$50/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 2:($50/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Polarization unset label @@ -398,7 +411,11 @@ file_in(mm) u 2:($35/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ file_in(mm) u 2:(-$41/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ file_in(mm) u 2:($41/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ file_in(mm) u 2:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 2:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 2:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 2:(-$43/sgn($5)) w l lc rgb nm1_color dt 1 title '$\gamma_p^{n=- 1}>0$' ,\ +file_in(mm) u 2:($43/sgn($5)) w l lc rgb nm1_color dt 2 title '$\gamma_p^{n=- 1}<0$' + + #Electron Damping Rates unset label @@ -413,10 +430,12 @@ set yrange [gam_min:gam_max] plot \ file_in(mm) u 2:(-$36/sgn($5)) w l lc rgb base_color dt 1 title '' ,\ file_in(mm) u 2:($36/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ -file_in(mm) u 2:(-$47/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ -file_in(mm) u 2:($47/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ -file_in(mm) u 2:(-$48/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 2:($48/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 2:(-$48/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ +file_in(mm) u 2:($48/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ +file_in(mm) u 2:(-$49/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ +file_in(mm) u 2:($49/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 2:(-$50/sgn($5)) w l lc rgb nm1_color dt 1 title '$\gamma_e^{n=- 1}>0$' ,\ +file_in(mm) u 2:($50/sgn($5)) w l lc rgb nm1_color dt 2 title '$\gamma_e^{n=- 1}>0$' #Polarization unset label diff --git a/plotter/example/example_kperp.plt b/plotter/example/example_kperp.plt index 06b4ad9..71ea504 100644 --- a/plotter/example/example_kperp.plt +++ b/plotter/example/example_kperp.plt @@ -43,6 +43,7 @@ dy=0.225 base_color='#000000' n0_color='#009e73' n1_color='#e68f00' +nm1_color='#00f86e' set log x set xrange [k_min:k_max] @@ -135,7 +136,9 @@ file_in(mm) u 1:($35/sgn($5)) w l lc rgb base_color dt 2 title '$\gamma_p<0$' ,\ file_in(mm) u 1:(-$41/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ file_in(mm) u 1:($41/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ file_in(mm) u 1:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 1:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 1:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 1:(-$43/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 1:($43/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Electron Damping Rates unset label @@ -151,10 +154,12 @@ set label '$|\gamma_e|/|\omega_{\textrm{r}}|$' at graph -0.3,0.15 rotate by 90 plot \ file_in(mm) u 1:(-$36/sgn($5)) w l lc rgb base_color dt 1 title '$\gamma_e>0$' ,\ file_in(mm) u 1:($36/sgn($5)) w l lc rgb base_color dt 2 title '$\gamma_e<0$' ,\ -file_in(mm) u 1:(-$47/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ -file_in(mm) u 1:($47/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ -file_in(mm) u 1:(-$48/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 1:($48/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 1:(-$48/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ +file_in(mm) u 1:($48/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ +file_in(mm) u 1:(-$49/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ +file_in(mm) u 1:($49/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 1:(-$50/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 1:($50/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Polarization unset label @@ -224,7 +229,9 @@ file_in(mm) u 1:($35/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ file_in(mm) u 1:(-$41/sgn($5)) w l lc rgb n0_color dt 1 title '$\gamma_p^{n=0}>0$' ,\ file_in(mm) u 1:($41/sgn($5)) w l lc rgb n0_color dt 2 title '$\gamma_p^{n=0}<0$' ,\ file_in(mm) u 1:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 1:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 1:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 1:(-$43/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 1:($43/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Electron Damping Rates unset label @@ -239,10 +246,12 @@ set yrange [gam_min:gam_max] plot \ file_in(mm) u 1:(-$36/sgn($5)) w l lc rgb base_color dt 1 title '' ,\ file_in(mm) u 1:($36/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ -file_in(mm) u 1:(-$47/sgn($5)) w l lc rgb n0_color dt 1 title '$\gamma_e^{n=0}>0$' ,\ -file_in(mm) u 1:($47/sgn($5)) w l lc rgb n0_color dt 2 title '$\gamma_e^{n=0}<0$' ,\ -file_in(mm) u 1:(-$48/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 1:($48/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 1:(-$48/sgn($5)) w l lc rgb n0_color dt 1 title '$\gamma_e^{n=0}>0$' ,\ +file_in(mm) u 1:($48/sgn($5)) w l lc rgb n0_color dt 2 title '$\gamma_e^{n=0}<0$' ,\ +file_in(mm) u 1:(-$49/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ +file_in(mm) u 1:($49/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 1:(-$50/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 1:($50/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Polarization unset label @@ -310,8 +319,10 @@ file_in(mm) u 1:(-$35/sgn($5)) w l lc rgb base_color dt 1 title '' ,\ file_in(mm) u 1:($35/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ file_in(mm) u 1:(-$41/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ file_in(mm) u 1:($41/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ -file_in(mm) u 1:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '$\gamma_p^{n=\pm 1}>0$' ,\ -file_in(mm) u 1:($42/sgn($5)) w l lc rgb n1_color dt 2 title '$\gamma_p^{n=\pm 1}<0$' +file_in(mm) u 1:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '$\gamma_p^{n=+ 1}>0$' ,\ +file_in(mm) u 1:($42/sgn($5)) w l lc rgb n1_color dt 2 title '$\gamma_p^{n=+ 1}<0$' ,\ +file_in(mm) u 1:(-$43/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 1:($43/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Electron Damping Rates unset label @@ -326,10 +337,12 @@ set yrange [gam_min:gam_max] plot \ file_in(mm) u 1:(-$36/sgn($5)) w l lc rgb base_color dt 1 title '' ,\ file_in(mm) u 1:($36/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ -file_in(mm) u 1:(-$47/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ -file_in(mm) u 1:($47/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ -file_in(mm) u 1:(-$48/sgn($5)) w l lc rgb n1_color dt 1 title '$\gamma_e^{n=\pm 1}>0$' ,\ -file_in(mm) u 1:($48/sgn($5)) w l lc rgb n1_color dt 2 title '$\gamma_e^{n=\pm 1}<0$' +file_in(mm) u 1:(-$48/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ +file_in(mm) u 1:($48/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ +file_in(mm) u 1:(-$49/sgn($5)) w l lc rgb n1_color dt 1 title '$\gamma_e^{n=+ 1}>0$' ,\ +file_in(mm) u 1:($49/sgn($5)) w l lc rgb n1_color dt 2 title '$\gamma_e^{n=+ 1}<0$' ,\ +file_in(mm) u 1:(-$50/sgn($5)) w l lc rgb nm1_color dt 1 title '' ,\ +file_in(mm) u 1:($50/sgn($5)) w l lc rgb nm1_color dt 2 title '' #Polarization unset label @@ -398,7 +411,9 @@ file_in(mm) u 1:($35/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ file_in(mm) u 1:(-$41/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ file_in(mm) u 1:($41/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ file_in(mm) u 1:(-$42/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 1:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 1:($42/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 2:(-$43/sgn($5)) w l lc rgb nm1_color dt 1 title '$\gamma_p^{n=- 1}>0$' ,\ +file_in(mm) u 2:($43/sgn($5)) w l lc rgb nm1_color dt 2 title '$\gamma_p^{n=- 1}>0$' #Electron Damping Rates unset label @@ -413,10 +428,12 @@ set yrange [gam_min:gam_max] plot \ file_in(mm) u 1:(-$36/sgn($5)) w l lc rgb base_color dt 1 title '' ,\ file_in(mm) u 1:($36/sgn($5)) w l lc rgb base_color dt 2 title '' ,\ -file_in(mm) u 1:(-$47/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ -file_in(mm) u 1:($47/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ -file_in(mm) u 1:(-$48/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ -file_in(mm) u 1:($48/sgn($5)) w l lc rgb n1_color dt 2 title '' +file_in(mm) u 1:(-$4/sgn($5)) w l lc rgb n0_color dt 1 title '' ,\ +file_in(mm) u 1:($48/sgn($5)) w l lc rgb n0_color dt 2 title '' ,\ +file_in(mm) u 1:(-$49/sgn($5)) w l lc rgb n1_color dt 1 title '' ,\ +file_in(mm) u 1:($49/sgn($5)) w l lc rgb n1_color dt 2 title '' ,\ +file_in(mm) u 1:(-$50/sgn($5)) w l lc rgb nm1_color dt 1 title '$\gamma_e^{n=- 1}>0$' ,\ +file_in(mm) u 1:($50/sgn($5)) w l lc rgb nm1_color dt 2 title '$\gamma_e^{n=- 1}>0$' #Polarization unset label diff --git a/src/disprels.f90 b/src/disprels.f90 index 99f2cc0..1bd7f01 100644 --- a/src/disprels.f90 +++ b/src/disprels.f90 @@ -748,15 +748,12 @@ subroutine om_scan(is) !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. - + !>>>GGH: 1/18/23 - real, dimension(1:6,1:nspec) :: Ps_split_new + real, dimension(1:7,1:nspec) :: Ps_split !!Power into/out of species/componets. !!Corrected LD/TTD calculation (GGH). + !!Split \pm 1 CD (KGK) real :: Ew !!Wave energy. @@ -892,15 +889,13 @@ subroutine om_scan(is) if ((scan(is)%eigen_s).or.((scan(is)%heat_s))) then val=abs(disp(omega)) !>>>GGH: 1/18/23 - call calc_eigen(omega,ef,bf,Us,ns,Ps,Ps_split,& - Ps_split_new,scan(is)%eigen_s,scan(is)%heat_s) + call calc_eigen(omega,ef,bf,Us,ns,Ps,& + Ps_split,scan(is)%eigen_s,scan(is)%heat_s) !<<>>GGH: 1/18/23 - Ps_split_new=-0.1 !<<>>GGH: 1/18/23 - real, dimension(1:6,1:nspec) :: Ps_split_new + real, dimension(1:7,1:nspec) :: Ps_split !!Power into/out of species/componets. !!Corrected LD/TTD calculation (GGH). + !!Split \pm 1 CD (KGK) real :: Ew !!Wave energy. @@ -1356,14 +1348,12 @@ subroutine om_double_scan if ((scan(2)%eigen_s).or.((scan(2)%heat_s))) then val=abs(disp(omega)) !>>>GGH: 1/18/23 - call calc_eigen(omega,ef,bf,Us,ns,Ps,Ps_split,Ps_split_new,scan(2)%eigen_s,scan(2)%heat_s) + call calc_eigen(omega,ef,bf,Us,ns,Ps,Ps_split,scan(2)%eigen_s,scan(2)%heat_s) !<<>>GGH: 1/18/23 - Ps_split_new=1.0 - !<<>>GGH: 1/18/23 - real, dimension(1:6,1:nspec) :: Ps_split_new + real, dimension(1:7,1:nspec) :: Ps_split !!Power into/out of species/componets. !!Corrected LD/TTD calculation (GGH). + !!Split \pm 1 CD (KGK) real :: Ew !!Wave energy. @@ -2087,7 +2063,7 @@ subroutine om_radial(omlast,params,out_unit,fmt,out_type,ir,mod_write) !Calculate eigenfunctions and heating rates if ((radial_heating).or.(radial_eigen)) then val=abs(disp(omega)) - call calc_eigen(omega,ef,bf,Us,ns,Ps,Ps_split,Ps_split_new,radial_eigen,radial_heating) + call calc_eigen(omega,ef,bf,Us,ns,Ps,Ps_split,radial_eigen,radial_heating) endif !Output results @@ -2126,7 +2102,7 @@ end subroutine om_radial !------------------------------------------------------------------------------ !-=-=-=-=-=- subroutine calc_eigen(omega,electric,magnetic,vmean,ns,Ps,& - Ps_split,Ps_split_new,eigen_L,heat_L) + Ps_split,eigen_L,heat_L) !! Calculates the electric and magnetic fields as well as species !! velocities and density fluctuations for identified wave !! as well as the power emission or absorption. @@ -2165,12 +2141,10 @@ subroutine calc_eigen(omega,electric,magnetic,vmean,ns,Ps,& real, dimension(1:nspec), intent(out) :: Ps !!Power into/out of species/components. - real, dimension(1:4,1:nspec), intent(out) :: Ps_split - !!Power into/out of species - - real, dimension(1:6,1:nspec), intent(out) :: Ps_split_new + real, dimension(1:7,1:nspec), intent(out) :: Ps_split !!<<>>GGH: 1/18/23 - write(fmt,'(a,i0,a)')'(6es15.6e3,',13*nspec,'es15.6e3)' + write(fmt,'(a,i0,a)')'(6es15.6e3,',14*nspec,'es15.6e3)' else write(fmt,'(a,i0,a)')'(6es15.6e3,',7*nspec,'es15.6e3)' endif @@ -3775,52 +3773,109 @@ complex function disp(om) !-=-=-=-=-=-=-=-=-=-= !n=\pm 1 susceptibility if ((low_n).and.(n==1)) then + + !\chi_xx susc_low(is,1,1,1) = norm(1)*& - n*n*jn(n,is)*( 2.*(disp_alph-1.) & + n*n*jn(n,is)*( (disp_alph-1.) & + + (sqrt(alphp * disp_tau/ disp_mu )/ kpar )*(& + disp_alph* (om - Vdrifts) * (zz(n)) & + + (zz(n))*(1.-disp_alph)*(n*disp_mu/(disp_Q )))) + + susc_low(is,1,1,-1) = norm(1)*& + n*n*jn(n,is)*( (disp_alph-1.) & + (sqrt(alphp * disp_tau/ disp_mu )/ kpar )*(& - disp_alph* (om - Vdrifts) * (zz(n)+zz(-n)) & - + (zz(n)-zz(-n))*(1.-disp_alph)*(n*disp_mu/(disp_Q )))) + disp_alph* (om - Vdrifts) * (zz(-n)) & + + (-zz(-n))*(1.-disp_alph)*(n*disp_mu/(disp_Q )))) + !\chi_xy susc_low(is,1,2,1) = norm(4)*& n*(jn(n,is)-jpn(n,is))*( & - (zz(n)-zz(-n))*disp_alph*(om-Vdrifts) & - +(zz(n)+zz(-n))*(1.-disp_alph)*n*disp_mu/(disp_Q )) + (zz(n))*disp_alph*(om-Vdrifts) & + +(zz(n))*(1.-disp_alph)*n*disp_mu/(disp_Q )) + + susc_low(is,1,2,-1) = norm(4)*& + n*(jn(n,is)-jpn(n,is))*( & + (-zz(-n))*disp_alph*(om-Vdrifts) & + +(zz(-n))*(1.-disp_alph)*n*disp_mu/(disp_Q )) + !\chi_yx susc_low(is,2,1,1) = -susc_low(is,1,2,1) + susc_low(is,2,1,-1) = -susc_low(is,1,2,-1) + + !\chi_xz susc_low(is,1,3,1) = norm(5)*& - n*jn(n,is)*(2.*n*(1.-disp_alph)*disp_mu/(disp_Q) + & + n*jn(n,is)*(n*(1.-disp_alph)*disp_mu/(disp_Q) + & (sqrt(alphp*disp_tau/disp_mu)/kpar)*(& - (zz(n) -zz(-n))*(disp_alph*om*(om-Vdrifts)-& + (zz(n))*(disp_alph*om*(om-Vdrifts)-& (n*n*disp_mu**2./(disp_Q**2.))*(1.-disp_alph)) + & - (zz(n) +zz(-n))*(disp_mu*n*(om-2.*om*disp_alph+Vdrifts*disp_alph)/disp_Q) )) + (zz(n))*(disp_mu*n*(om-2.*om*disp_alph+Vdrifts*disp_alph)/disp_Q) )) + susc_low(is,1,3,-1) = norm(5)*& + n*jn(n,is)*(n*(1.-disp_alph)*disp_mu/(disp_Q) + & + (sqrt(alphp*disp_tau/disp_mu)/kpar)*(& + (-zz(-n))*(disp_alph*om*(om-Vdrifts)-& + (n*n*disp_mu**2./(disp_Q**2.))*(1.-disp_alph)) + & + (zz(-n))*(disp_mu*n*(om-2.*om*disp_alph+Vdrifts*disp_alph)/disp_Q) )) + + !\chi_zx susc_low(is,3,1,1) = susc_low(is,1,3,1) + susc_low(is,3,1,-1) = susc_low(is,1,3,-1) + !\chi_yy susc_low(is,2,2,1) = norm(2)*& (n*n*jn(n,is) + 2.*lambdas*lambdas*(jn(n,is)-jpn(n,is)))*& - (2.*(disp_alph-1.) + (sqrt(alphp*disp_tau/disp_mu)/kpar)*(& - disp_alph*(om - Vdrifts)*(zz(n)+zz(-n)) + (n*disp_mu/(disp_Q))*& - (1.-disp_alph)*(zz(n)-zz(-n)))) + ((disp_alph-1.) + (sqrt(alphp*disp_tau/disp_mu)/kpar)*(& + disp_alph*(om - Vdrifts)*(zz(n)) + (n*disp_mu/(disp_Q))*& + (1.-disp_alph)*(zz(n)))) + + susc_low(is,2,2,-1) = norm(2)*& + (n*n*jn(n,is) + 2.*lambdas*lambdas*(jn(n,is)-jpn(n,is)))*& + ((disp_alph-1.) + (sqrt(alphp*disp_tau/disp_mu)/kpar)*(& + disp_alph*(om - Vdrifts)*(zz(-n)) + (n*disp_mu/(disp_Q))*& + (1.-disp_alph)*(-zz(-n)))) + !\chi_yz susc_low(is,2,3,1) = norm(6)*& - (jn(n,is)-jpn(n,is))*(2.*(om*disp_alph-Vdrifts) & + (jn(n,is)-jpn(n,is))*((om*disp_alph-Vdrifts) & + +(sqrt(alphp*disp_tau/disp_mu)/kpar)*(& + (zz(n))*(disp_alph*om*(om-Vdrifts) -& + (n*n*disp_mu**2./(disp_Q**2.))*(1.-disp_alph)) & + +(zz(n))*(n*disp_mu*(om-disp_alph*(2.*om-Vdrifts))/disp_Q) )) + + susc_low(is,2,3,-1) = norm(6)*& + (jn(n,is)-jpn(n,is))*((om*disp_alph-Vdrifts) & +(sqrt(alphp*disp_tau/disp_mu)/kpar)*(& - (zz(n)+zz(-n))*(disp_alph*om*(om-Vdrifts) -& + (zz(-n))*(disp_alph*om*(om-Vdrifts) -& (n*n*disp_mu**2./(disp_Q**2.))*(1.-disp_alph)) & - +(zz(n)-zz(-n))*(n*disp_mu*(om-disp_alph*(2.*om-Vdrifts))/disp_Q) )) + +(-zz(-n))*(n*disp_mu*(om-disp_alph*(2.*om-Vdrifts))/disp_Q) )) - susc_low(is,3,2,1) = -susc_low(is,2,3,1) + !\chi_zy + susc_low(is,3,2,1) = -susc_low(is,2,3,1) + + susc_low(is,3,2,-1) = -susc_low(is,2,3,-1) + !\chi_zz susc_low(is,3,3,1) = norm(3)*& jn(n,is)*(2.*om* (om*disp_alph - Vdrifts) - & 2.*n*n*disp_mu**2.* (1.-disp_alph)/(disp_Q**2.) & +(sqrt(alphp*disp_tau/disp_mu)/kpar) * (& - (zz(n)+zz(-n))*(om*om*disp_alph * (om - Vdrifts) & + (zz(n))*(om*om*disp_alph * (om - Vdrifts) & + n*n*disp_mu**2.*(om*(3.*disp_alph-2.) -disp_alph*Vdrifts)/disp_Q**2.) & - +(n * disp_mu/disp_Q )*(zz(n)-zz(-n))*& !changed sign + +(n * disp_mu/disp_Q )*(zz(n))*& !changed sign (om*(om+2.*disp_alph *Vdrifts -3.*disp_alph*om)+& n*n*disp_mu**2.*(1.-disp_alph)/( disp_Q**2.)) ) )!? + + susc_low(is,3,3,-1) = norm(3)*& + jn(n,is)*(2.*om* (om*disp_alph - Vdrifts) - & + 2.*n*n*disp_mu**2.* (1.-disp_alph)/(disp_Q**2.) & + +(sqrt(alphp*disp_tau/disp_mu)/kpar) * (& + (zz(-n))*(om*om*disp_alph * (om - Vdrifts) & + + n*n*disp_mu**2.*(om*(3.*disp_alph-2.) -disp_alph*Vdrifts)/disp_Q**2.) & + +(n * disp_mu/disp_Q )*(-zz(-n))*& !changed sign + (om*(om+2.*disp_alph *Vdrifts -3.*disp_alph*om)+& + n*n*disp_mu**2.*(1.-disp_alph)/( disp_Q**2.)) ) )!? + endif enddo diff --git a/src/fpc.f90 b/src/fpc.f90 index fd7e7f6..86b7afe 100644 --- a/src/fpc.f90 +++ b/src/fpc.f90 @@ -106,10 +106,11 @@ subroutine compute_fpc_cart(wrootindex) real, dimension(1:nspec) :: Ps !! Power into/out of species - real, dimension(1:4, 1:nspec) :: Ps_split + !real, dimension(1:4, 1:nspec) :: Ps_split !! Power into/out of species (Tensor that holds different channels (TTD, LD, CD)) + !! Deprecated - real, dimension(1:6, 1:nspec) :: Ps_split_new + real, dimension(1:7, 1:nspec) :: Ps_split !! Power into/out of species updated by Greg G Howes to include off diagnal components real :: Ew @@ -293,7 +294,7 @@ subroutine compute_fpc_cart(wrootindex) iflag = 0 omega = rtsec(disp, om1, om2, tol, iflag) - call calc_eigen(omega, ef, bf, Us, ns, Ps, Ps_split, Ps_split_new, .true., .true.) + call calc_eigen(omega, ef, bf, Us, ns, Ps, Ps_split, .true., .true.) if (ABS(aimag(omega)) .gt. 1./eeuler) then write (*, *) 'WARNING: damping term of omega is stronger than 1/eeuler' @@ -846,33 +847,36 @@ subroutine compute_fpc_cart(wrootindex) if (ABS(Ps(is)) .lt. 9.999E-99) then Ps(is) = 0. end if - if (ABS(Ps_split_new(1,is)) .lt. 9.999E-99) then - Ps_split_new(1,is) = 0. + if (ABS(Ps_split(1,is)) .lt. 9.999E-99) then + Ps_split(1,is) = 0. end if - if (ABS(Ps_split_new(2,is)) .lt. 9.999E-99) then - Ps_split_new(2,is) = 0. + if (ABS(Ps_split(2,is)) .lt. 9.999E-99) then + Ps_split(2,is) = 0. end if - if (ABS(Ps_split_new(3,is)) .lt. 9.999E-99) then - Ps_split_new(3,is) = 0. + if (ABS(Ps_split(3,is)) .lt. 9.999E-99) then + Ps_split(3,is) = 0. end if - if (ABS(Ps_split_new(4,is)) .lt. 9.999E-99) then - Ps_split_new(4,is) = 0. + if (ABS(Ps_split(4,is)) .lt. 9.999E-99) then + Ps_split(4,is) = 0. end if - if (ABS(Ps_split_new(5,is)) .lt. 9.999E-99) then - Ps_split_new(5,is) = 0. + if (ABS(Ps_split(5,is)) .lt. 9.999E-99) then + Ps_split(5,is) = 0. end if - if (ABS(Ps_split_new(6,is)) .lt. 9.999E-99) then - Ps_split_new(6,is) = 0. + if (ABS(Ps_split(6,is)) .lt. 9.999E-99) then + Ps_split(6,is) = 0. + endif + if (ABS(Ps_split(7,is)) .lt. 9.999E-99) then + Ps_split(7,is) = 0. end if end do !Write format (consistent with usual PLUME output) - write (fmt, '(a,i0,a)') '(6es15.6,12es15.6,', 15*nspec, 'es15.6)' + write (fmt, '(a,i0,a)') '(6es15.6,12es15.6,', 16*nspec, 'es15.6)' write (unit_s + 5, fmt) & kperp, kpar, betap, vtp, & omega, & bf(1:3), ef(1:3), Us(1:3, 1:nspec), ns(1:nspec), & - Ps(1:nspec), Ps_split_new(1:6, 1:nspec) + Ps(1:nspec), Ps_split(1:7, 1:nspec) close (unit_s + 5) !Write Velocity Integrated Moments------------------------------------------- @@ -1014,10 +1018,11 @@ subroutine compute_fpc_gyro(wrootindex) real, dimension(1:nspec) :: Ps !! Power into/out of species - real, dimension(1:4, 1:nspec) :: Ps_split + !real, dimension(1:4, 1:nspec) :: Ps_split !! Power into/out of species + !! Deprecated - real, dimension(1:6, 1:nspec) :: Ps_split_new + real, dimension(1:7, 1:nspec) :: Ps_split !! Power into/out of species (includes all diagonal terms added by GGH in 2023) real :: Ew @@ -1145,7 +1150,7 @@ subroutine compute_fpc_gyro(wrootindex) iflag = 0 omega = rtsec(disp, om1, om2, tol, iflag) - call calc_eigen(omega, ef, bf, Us, ns, Ps, Ps_split, Ps_split_new, .true., .true.) + call calc_eigen(omega, ef, bf, Us, ns, Ps, Ps_split, .true., .true.) if (ABS(aimag(omega)) .gt. 1./eeuler) then write (*, *) 'WARNING: damping term of omega is stronger than 1/eeuler- the assumptions made to derive the FPC ' diff --git a/src/functions.f90 b/src/functions.f90 index f7ed4b5..e5653a6 100644 --- a/src/functions.f90 +++ b/src/functions.f90 @@ -86,7 +86,7 @@ subroutine read_in_params !Allocate the susceptibility tensor to have nspec indicies allocate (susc(1:nspec, 3, 3)) if (low_n) & - allocate (susc_low(1:nspec, 3, 3, 0:1)) + allocate (susc_low(1:nspec, 3, 3, -1:1)) !initialized quasineutrality and current conservation check. sum_nq = 0.; sum_nqv = 0. diff --git a/src/vars.f90 b/src/vars.f90 index c979841..b778480 100644 --- a/src/vars.f90 +++ b/src/vars.f90 @@ -167,7 +167,7 @@ module vars complex, dimension(:, :, :, :), allocatable:: susc_low !! low-n components of the susceptibility tensor. - !! (1:nspec,1:3,1:3,0:1) with the 3x3 subarray arranged as susc. + !! (1:nspec,1:3,1:3,-1:1) with the 3x3 subarray arranged as susc. !! The final index contains the n=0 and n=\pm 1 contributions. logical :: low_n = .true.