Changeset 40491

Timestamp:

Jun 27, 2018, 4:11:56 PM (8 years ago)

Author:

eugene

Message:

add new psf model type for HSC

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/objects/Makefile.am (modified) (2 diffs)
• src/objects/models/fwhm.sh (modified) (7 diffs)
• src/objects/models/hscv1.fwhm.h (copied) (copied from branches/czw_branch/20170908/psModules/src/objects/models/hscv1.fwhm.h )
• src/objects/models/pmModel_HSC_V1.c (copied) (copied from branches/czw_branch/20170908/psModules/src/objects/models/pmModel_HSC_V1.c )
• src/objects/models/pmModel_HSC_V1.h (copied) (copied from branches/czw_branch/20170908/psModules/src/objects/models/pmModel_HSC_V1.h )
• src/objects/pmModelClass.c (modified) (2 diffs)
• src/objects/pmSourceOutputs.c (modified) (1 diff)

trunk/psModules/src/objects/Makefile.am

@@ -92 +92 @@
         models/pmModel_EXP.c \
         models/pmModel_DEV.c \
-        models/pmModel_TRAIL.c
+        models/pmModel_TRAIL.c \
+        models/pmModel_HSC_V1.c
 
 pkginclude_HEADERS = \
@@ -142 +143 @@
         models/pmModel_EXP.h \
         models/pmModel_DEV.h \
-        models/pmModel_TRAIL.h
+        models/pmModel_TRAIL.h \
+        models/pmModel_HSC_V1.h
 
 CLEANFILES = *~

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

@@ -161 +161 @@
 end
 
-# qgauss is like ps1_v1 with z^2.25
-macro find.fwhm.qgauss
+# hsc
+macro find.fwhm.hscv1
   if ($0 != 2)
-    echo "USAGE: find.qgauss (K)"
+    echo "USAGE: find.fwhm (K)"
     break
   end
@@ -180 +180 @@
   create z 0 50 $dz
 
-  # set f = 1.0 / (1.0 + $K*z + z^2.25)
-  # f = $K*z + z^2.25 - 1.0, find where f == 0.0
-  
-  set f = $K*z + z^2.25 - 1.0
-  set dfdz = $K + 2.25*z^1.25
+  # set f = 1.0 / (1.0 + $K*z + z^1.667)
+  # f = $K*z + z^1.667 - 1.0, find where f == 0.0
+  
+  set f = $K*z + z^1.8 - 1.0
+  set dfdz = $K + 1.8*z^0.8
 
   #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
+  # these just need to bound the solution
+  $nZ0 = 0.5 / $dz
+  $nZ1 = 5.0 / $dz
 
   $Zg = 0.5*(z[$nZ0] + z[$nZ1])  
@@ -206 +207 @@
     $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
+  $Zhm = $Zg
+  $FWHM = 2*sqrt(2*$Zg)
+
+  fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm 
+end
+
+# qgauss is like ps1_v1 with z^2.25
+macro find.fwhm.qgauss
+  if ($0 != 2)
+    echo "USAGE: find.qgauss (K)"
+    break
+  end
+  
+  $K = $1
+
+  if ($K == 0.0)
+    $Zhm = 1.0
+    $FWHM = 2*sqrt(2)
+    fprintf "%4.2f, // %4.1f, %4.2f" $FWHM $K $Zhm 
+    return
+  end
+
+  $dz = 0.01
+  create z 0 50 $dz
+
+  # set f = 1.0 / (1.0 + $K*z + z^2.25)
+  # f = $K*z + z^2.25 - 1.0, find where f == 0.0
+  
+  set f = $K*z + z^2.25 - 1.0
+  set dfdz = $K + 2.25*z^1.25
+
+  #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
+    # 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 $Zg $Fg $dZ $dFdz_g
 
@@ -223 +286 @@
   if ($0 != 4)
     echo "USAGE: fwhm.trend (model) (struct) (output)"
-    echo "  model: pgauss, rgauss, ps1v1, qgauss"
+    echo "  model: pgauss, rgauss, ps1v1, qgauss, hscv1"
     break
   end
@@ -229 +292 @@
   delete fwhm_v k_v
 
-  $SAVE = 0
+  $SAVE = 1
   if ($SAVE) exec rm -f $3
   if ($SAVE) output $3
@@ -240 +303 @@
   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[] = \{"
+#  if ($SAVE) echo "static float $2[] = \{"
 
   for k $minK $maxK $delK -incl
@@ -247 +310 @@
     concat $FWHM fwhm_v
   end
-  if ($SAVE) echo "\}@"
+#  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
-
+#  lim k_v fwhm_v; clear; box; plot k_v fwhm_v -pt 10 -sz 1.0
+end
+

trunk/psModules/src/objects/pmModelClass.c

@@ -47 +47 @@
 # include "models/pmModel_QGAUSS.h"
 # include "models/pmModel_PS1_V1.h"
+# include "models/pmModel_HSC_V1.h"
 # include "models/pmModel_RGAUSS.h"
 # include "models/pmModel_SERSIC.h"
@@ -63 +64 @@
     {"PS_MODEL_DEV",          7, 1, (pmModelFunc)pmModelFunc_DEV,     (pmModelFlux)pmModelFlux_DEV,     (pmModelRadius)pmModelRadius_DEV,     (pmModelSetFWHM)pmModelSetFWHM_DEV,     (pmModelLimits)pmModelLimits_DEV,     (pmModelGuessFunc)pmModelGuess_DEV,    (pmModelFromPSFFunc)pmModelFromPSF_DEV,    (pmModelParamsFromPSF)pmModelParamsFromPSF_DEV,    (pmModelFitStatusFunc)pmModelFitStatus_DEV,    (pmModelSetLimitsFunc)pmModelSetLimits_DEV    },
     {"PS_MODEL_TRAIL",        7, 0, (pmModelFunc)pmModelFunc_TRAIL,   (pmModelFlux)pmModelFlux_TRAIL,   (pmModelRadius)pmModelRadius_TRAIL,   (pmModelSetFWHM)pmModelSetFWHM_TRAIL,   (pmModelLimits)pmModelLimits_TRAIL,   (pmModelGuessFunc)pmModelGuess_TRAIL,  (pmModelFromPSFFunc)pmModelFromPSF_TRAIL,  (pmModelParamsFromPSF)pmModelParamsFromPSF_TRAIL,  (pmModelFitStatusFunc)pmModelFitStatus_TRAIL,  (pmModelSetLimitsFunc)pmModelSetLimits_TRAIL  },
+    {"PS_MODEL_HSC_V1",       8, 0, (pmModelFunc)pmModelFunc_HSC_V1,  (pmModelFlux)pmModelFlux_HSC_V1,  (pmModelRadius)pmModelRadius_HSC_V1,  (pmModelSetFWHM)pmModelSetFWHM_HSC_V1,  (pmModelLimits)pmModelLimits_HSC_V1,  (pmModelGuessFunc)pmModelGuess_HSC_V1, (pmModelFromPSFFunc)pmModelFromPSF_HSC_V1, (pmModelParamsFromPSF)pmModelParamsFromPSF_HSC_V1, (pmModelFitStatusFunc)pmModelFitStatus_HSC_V1, (pmModelSetLimitsFunc)pmModelSetLimits_HSC_V1 },
 };
 

trunk/psModules/src/objects/pmSourceOutputs.c

@@ -135 +135 @@
                 outputs->psfCore = PAR[PM_PAR_7];
             }
+            if (model->type == pmModelClassGetType ("PS_MODEL_HSC_V1")) {
+                outputs->psfCore = PAR[PM_PAR_7];
+            }
 
             outputs->psfMajorFWHM = model->class->modelSetFWHM(model->params, axes.major);