From 2e817f95d42c773ff0181db4c34370ebac7fb6e0 Mon Sep 17 00:00:00 2001
From: William Henney <whenney@gmail.com>
Date: Wed, 17 Nov 2010 18:10:01 -0600
Subject: [PATCH] Added scaling to  cubedenstats.f90 to prevent FP overflow

---
 .gitignore       |  2 ++
 cubedenstats.f90 | 11 +++++++++++
 2 files changed, 13 insertions(+)

diff --git a/.gitignore b/.gitignore
index 0c7f1b9..6bf9afd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -17,3 +17,5 @@
 /phabrebin
 /rotatetest
 /turbstats
+/cubeextras
+/cubet
diff --git a/cubedenstats.f90 b/cubedenstats.f90
index cc63d6a..b123394 100644
--- a/cubedenstats.f90
+++ b/cubedenstats.f90
@@ -79,6 +79,8 @@ program cubedenstats
      dmean_tot = sum(d)/real(nx*ny*nz, dp)
      d2mean_tot = sum(d*d)/real(nx*ny*nz, dp)
 
+     d = d / dmean_tot          ! rescale density to prevent overflow
+
      ! version 2: weight by the ion fraction
      dmean_m = sum(d*wm)/sum(wm)
      dmean_n = sum(d*wn)/sum(wn)
@@ -92,6 +94,15 @@ program cubedenstats
      ! mean of density-cubed - only for ionized gas
      d3mean_i = sum(d*d*d*wi, mask=m)/sum(wi, mask=m)
 
+     ! put scaling back in
+     dmean_m = dmean_m * dmean_tot
+     dmean_n = dmean_n * dmean_tot
+     dmean_i = dmean_i * dmean_tot
+     d2mean_m = d2mean_m * dmean_tot * dmean_tot
+     d2mean_n = d2mean_n * dmean_tot * dmean_tot
+     d2mean_i = d2mean_i * dmean_tot * dmean_tot
+     d3mean_i = d3mean_i * dmean_tot * dmean_tot * dmean_tot
+
      if (mod(it,10)==0) print *, 'Done timestep: ', it
 
      write(1, '(i4.4,"'//TAB//'",9(es11.3,"'//TAB//'")))') it, &