Merge pull request #597 from MESAHub/update_sneut

replace sin/cos with multiangle formula in neutrino cooling (part II)
MESAHub · Dec 11, 2023 · 39298f7 · 39298f7
2 parents a563abf + 8abfe95
commit 39298f7
Show file tree

Hide file tree

Showing 6 changed files with 41 additions and 10,019 deletions.
diff --git a/.gitattributes b/.gitattributes
@@ -10,6 +10,7 @@
 *.mp4 filter=lfs diff=lfs merge=lfs -text
 data/star_data/zams_models/zams_z2m2_y28.data filter=lfs diff=lfs merge=lfs -text
 docs/source/assets/zams_z2m2_y28_patched.data filter=lfs diff=lfs merge=lfs -text
+neu/test/test_output filter=lfs diff=lfs merge=lfs -text
 
 # Try to get github to better understand what our filetypes are
 *.inc linguist-language=Fortran

diff --git a/gyre/private/.gitignore b/gyre/private/.gitignore
diff --git a/linters/__init__.py b/linters/__init__.py
diff --git a/neu/private/mod_neu.f90 b/neu/private/mod_neu.f90
@@ -694,17 +694,26 @@ subroutine phot_neu(sphot,sphotdt,sphotdd,sphotda,sphotdz, input)
 
    !..equation 3.7, compute the expensive trig functions only one time
          cos1 = cos(fac1*tau)
-         cos2 = cos(fac1*2.0d0*tau)
-         cos3 = cos(fac1*3.0d0*tau)
-         cos4 = cos(fac1*4.0d0*tau)
-         cos5 = cos(fac1*5.0d0*tau)
-         last = cos(fac2*tau)
-
          sin1 = sin(fac1*tau)
-         sin2 = sin(fac1*2.0d0*tau)
-         sin3 = sin(fac1*3.0d0*tau)
-         sin4 = sin(fac1*4.0d0*tau)
-         sin5 = sin(fac1*5.0d0*tau)
+
+         ! double, triple, etc. angle formulas
+         ! sin/cos (2 fac1 tau)
+         sin2 = 2.0d0 * sin1 * cos1
+         cos2 = 2.0d0 * cos1 * cos1 - 1.0d0
+
+         ! sin/cos (3 fac1 tau)
+         sin3 = sin1 * (3.0d0 - 4.0d0 * sin1 * sin1)
+         cos3 = cos1 * (4.0d0 * cos1 * cos1 - 3.0d0)
+
+         ! sin/cos (4 fac1 tau) -- use double angle on sin2/cos2
+         sin4 = 2.0d0 * sin2 * cos2
+         cos4 = 2.0d0 * cos2 * cos2 - 1.0d0
+
+         ! sin/cos (5 fac1 tau)
+         sin5 = sin1 * (5.0d0 - sin1 * sin1 * (20.0d0 - 16.0d0 * sin1 * sin1))
+         cos5 = cos1 * (cos1 * cos1 * (16.0d0 * cos1 * cos1 - 20.0d0) + 5.0d0)
+
+         last = cos(fac2*tau)
          xast = sin(fac2*tau)
 
          a0 = 0.5d0*c00  &
@@ -1021,16 +1030,25 @@ subroutine brem_neu_liquid_metal(sbrem,sbremdt,sbremdd,sbremda,sbremdz,t8, input
 
    !..compute the expensive trig functions of equation 5.21 only once
          cos1 = cos(u)
-         cos2 = cos(2.0d0*u)
-         cos3 = cos(3.0d0*u)
-         cos4 = cos(4.0d0*u)
-         cos5 = cos(5.0d0*u)
-
          sin1 = sin(u)
-         sin2 = sin(2.0d0*u)
-         sin3 = sin(3.0d0*u)
-         sin4 = sin(4.0d0*u)
-         sin5 = sin(5.0d0*u)
+
+         ! double, triple, etc. angle formulas
+         ! sin/cos (2 u)
+         sin2 = 2.0d0 * sin1 * cos1
+         cos2 = 2.0d0 * cos1 * cos1 - 1.0d0
+
+         ! sin/cos (3 u)
+         sin3 = sin1 * (3.0d0 - 4.0d0 * sin1 * sin1)
+         cos3 = cos1 * (4.0d0 * cos1 * cos1 - 3.0d0)
+
+         ! sin/cos (4 u) -- use double angle on sin2/cos2
+         sin4 = 2.0d0 * sin2 * cos2
+         cos4 = 2.0d0 * cos2 * cos2 - 1.0d0
+
+         ! sin/cos (5 u)
+         sin5 = sin1 * (5.0d0 - sin1 * sin1 * (20.0d0 - 16.0d0 * sin1 * sin1))
+         cos5 = cos1 * (cos1 * cos1 * (16.0d0 * cos1 * cos1 - 20.0d0) + 5.0d0)
+
 
    !..equation 5.21
          fb =  0.5d0 * 0.17946d0  + 0.00945d0*u + 0.34529d0    &