Changeset 38280 for trunk/psModules/src/objects/models/fwhm.sh

Timestamp:

May 15, 2015, 1:56:42 PM (11 years ago)

Author:

eugene

Message:

fix interpolation of Sigma,Core -> FWHM in PS1_V1 and QGAUSS; test code to prove fix

File:

• trunk/psModules/src/objects/models/fwhm.sh (modified) (8 diffs)

trunk/psModules/src/objects/models/fwhm.sh

@@ -110 +110 @@
     $Zhm = 1.0
     $FWHM = 2*sqrt(2)
-    echo $K : $Zhm : $FWHM
+    fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm 
     return
   end
@@ -121 +121 @@
   
   set f = $K*z + z^1.667 - 1.0
-  set dfdz = $K + z^0.667
+  set dfdz = $K + 1.677*z^0.667
 
   #lim -n 0 z f; clear; box; plot z f
   #lim -n 1 z dfdz; clear; box; plot z dfdz  
 
-  $nZ0 = 0
-  $nZ1 = 5 / $dz
-
-  $Zg = 0.5*(z[$nZ0] + z[$nZ1])  
-  $nZg = int($Zg / $dz)
-  $dZ = 1.0
-
-  for i 0 10
-    $Fg = f[$nZg]
-    $dFdz_g = dfdz[$nZg]
-
-    $dZ = $Fg / $dFdz_g
-
-    # echo $Zg $Fg $dZ $dFdz_g
-
-    $Zg -= $dZ
-    $Zg = max ($Zg , 0)
-    $nZg = int($Zg / $dz)
+  # these just need to bound the solution
+  $nZ0 = 0.5 / $dz
+  $nZ1 = 5.0 / $dz
+
+  $Zg = 0.5*(z[$nZ0] + z[$nZ1])  
+  $nZg = int($Zg / $dz)
+  $dZ = 1.0
+
+  for i 0 10
+    # interpolate between $nZg and $nZg + 1
+    $Fg = f[$nZg] + ($Zg - $nZg*$dz)*(f[$nZg+1] - f[$nZg])/$dz
+    $dFdz_g = dfdz[$nZg] + ($Zg - $nZg*$dz)*(dfdz[$nZg+1] - dfdz[$nZg])/$dz
+
+    # $Fg = f[$nZg]
+    # $dFdz_g = dfdz[$nZg]
+
+    $dZ = $Fg / $dFdz_g
+
+    # echo -no-return $Zg $Fg $dZ $dFdz_g
+
+    $Zg -= $dZ
+
+    $Zg = max ($Zg , 0)
+    $Zg = min ($Zg , z[-2])
+    $nZg = int($Zg / $dz)
+    
+    # echo ": $Zg $nZg"
   end 
   # echo $Zg $Fg $dZ $dFdz_g
@@ -149 +158 @@
   $FWHM = 2*sqrt(2*$Zg)
 
-  echo $K : $Zhm : $FWHM
+  fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm 
 end
 
@@ -164 +173 @@
     $Zhm = 1.0
     $FWHM = 2*sqrt(2)
-    echo $K : $Zhm : $FWHM
+    fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm 
     return
   end
@@ -175 +184 @@
   
   set f = $K*z + z^2.25 - 1.0
-  set dfdz = $K + z^1.25
+  set dfdz = $K + 2.25*z^1.25
 
   #lim -n 0 z f; clear; box; plot z f
@@ -188 +197 @@
 
   for i 0 10
-    $Fg = f[$nZg]
-    $dFdz_g = dfdz[$nZg]
+    # interpolate between $nZg and $nZg + 1
+    $Fg = f[$nZg] + ($Zg - $nZg*$dz)*(f[$nZg+1] - f[$nZg])/$dz
+    $dFdz_g = dfdz[$nZg] + ($Zg - $nZg*$dz)*(dfdz[$nZg+1] - dfdz[$nZg])/$dz
+
+    # $Fg = f[$nZg]
+    # $dFdz_g = dfdz[$nZg]
 
     $dZ = $Fg / $dFdz_g
@@ -197 +210 @@
     $Zg -= $dZ
     $Zg = max ($Zg , 0)
+    $Zg = min ($Zg , z[-2])
     $nZg = int($Zg / $dz)
   end 
@@ -203 +217 @@
   $FWHM = 2*sqrt(2*$Zg)
 
-  echo $K : $Zhm : $FWHM
+  fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm 
 end
 
 macro fwhm.trend
+  if ($0 != 4)
+    echo "USAGE: fwhm.trend (model) (struct) (output)"
+    echo "  model: pgauss, rgauss, ps1v1, qgauss"
+    break
+  end
 
   delete fwhm_v k_v
 
-  for k 0 20 1.0
-    find.fwhm.qgauss $k
+  $SAVE = 0
+  if ($SAVE) exec rm -f $3
+  if ($SAVE) output $3
+
+  $minK = -1
+  $maxK = 20
+  $delK = 0.2
+
+  if ($SAVE) echo "# define FWHM_BIN $delK"
+  if ($SAVE) echo "# define MIN_FWHM_BIN $minK"
+  if ($SAVE) echo "# define N_FWHM_BIN {int(($maxK - $minK) / $delK) + 1}"
+  if ($SAVE) echo "static float $2[] = \{"
+
+  for k $minK $maxK $delK -incl
+    find.fwhm.$1 $k
     concat $k k_v
     concat $FWHM fwhm_v
   end
-end
-
+  if ($SAVE) echo "\}@"
+  if ($SAVE) output stdout
+
+  lim k_v fwhm_v; clear; box; plot k_v fwhm_v -pt 10 -sz 1.0
+end
+