Index: branches/czw_branch/20170908/ippconfig/dvo.photcodes =================================================================== --- branches/czw_branch/20170908/ippconfig/dvo.photcodes (revision 40129) +++ branches/czw_branch/20170908/ippconfig/dvo.photcodes (revision 40179) @@ -1663,3 +1663,115 @@ 25110 HSC.NB921.110 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 25111 HSC.NB921.111 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 - + 26000 HSC.NB816.00 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26001 HSC.NB816.01 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26002 HSC.NB816.02 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26003 HSC.NB816.03 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26004 HSC.NB816.04 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26005 HSC.NB816.05 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26006 HSC.NB816.06 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26007 HSC.NB816.07 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26008 HSC.NB816.08 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26009 HSC.NB816.09 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26010 HSC.NB816.10 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26011 HSC.NB816.11 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26012 HSC.NB816.12 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26013 HSC.NB816.13 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26014 HSC.NB816.14 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26015 HSC.NB816.15 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26016 HSC.NB816.16 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26017 HSC.NB816.17 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26018 HSC.NB816.18 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26019 HSC.NB816.19 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26020 HSC.NB816.20 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26021 HSC.NB816.21 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26022 HSC.NB816.22 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26023 HSC.NB816.23 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26024 HSC.NB816.24 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26025 HSC.NB816.25 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26026 HSC.NB816.26 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26027 HSC.NB816.27 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26028 HSC.NB816.28 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26029 HSC.NB816.29 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26030 HSC.NB816.30 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26031 HSC.NB816.31 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26032 HSC.NB816.32 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26033 HSC.NB816.33 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26034 HSC.NB816.34 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26035 HSC.NB816.35 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26036 HSC.NB816.36 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26037 HSC.NB816.37 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26038 HSC.NB816.38 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26039 HSC.NB816.39 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26040 HSC.NB816.40 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26041 HSC.NB816.41 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26042 HSC.NB816.42 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26043 HSC.NB816.43 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26044 HSC.NB816.44 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26045 HSC.NB816.45 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26046 HSC.NB816.46 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26047 HSC.NB816.47 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26048 HSC.NB816.48 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26049 HSC.NB816.49 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26050 HSC.NB816.50 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26051 HSC.NB816.51 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26052 HSC.NB816.52 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26053 HSC.NB816.53 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26054 HSC.NB816.54 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26055 HSC.NB816.55 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26056 HSC.NB816.56 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26057 HSC.NB816.57 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26058 HSC.NB816.58 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26059 HSC.NB816.59 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26060 HSC.NB816.60 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26061 HSC.NB816.61 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26062 HSC.NB816.62 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26063 HSC.NB816.63 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26064 HSC.NB816.64 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26065 HSC.NB816.65 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26066 HSC.NB816.66 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26067 HSC.NB816.67 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26068 HSC.NB816.68 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26069 HSC.NB816.69 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26070 HSC.NB816.70 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26071 HSC.NB816.71 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26072 HSC.NB816.72 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26073 HSC.NB816.73 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26074 HSC.NB816.74 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26075 HSC.NB816.75 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26076 HSC.NB816.76 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26077 HSC.NB816.77 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26078 HSC.NB816.78 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26079 HSC.NB816.79 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26080 HSC.NB816.80 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26081 HSC.NB816.81 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26082 HSC.NB816.82 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26083 HSC.NB816.83 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26084 HSC.NB816.84 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26085 HSC.NB816.85 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26086 HSC.NB816.86 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26087 HSC.NB816.87 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26088 HSC.NB816.88 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26089 HSC.NB816.89 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26090 HSC.NB816.90 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26091 HSC.NB816.91 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26092 HSC.NB816.92 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26093 HSC.NB816.93 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26094 HSC.NB816.94 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26095 HSC.NB816.95 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26096 HSC.NB816.96 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26097 HSC.NB816.97 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26098 HSC.NB816.98 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26099 HSC.NB816.99 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26100 HSC.NB816.100 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26101 HSC.NB816.101 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26102 HSC.NB816.102 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26103 HSC.NB816.103 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26104 HSC.NB816.104 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26105 HSC.NB816.105 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26106 HSC.NB816.106 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26107 HSC.NB816.107 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26108 HSC.NB816.108 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26109 HSC.NB816.109 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26110 HSC.NB816.110 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + 26111 HSC.NB816.111 dep 25.000 0.000 0.000 - - 0.0000 0 - 0.000 0.000 0.000 0.000 0x0000 0x0000 0x0000 0x0000 + Index: branches/czw_branch/20170908/ippconfig/hsc/Makefile.am =================================================================== --- branches/czw_branch/20170908/ippconfig/hsc/Makefile.am (revision 40129) +++ branches/czw_branch/20170908/ippconfig/hsc/Makefile.am (revision 40179) @@ -11,4 +11,5 @@ ppMerge.config \ psastro.config \ + psphot.config \ ppStats.config \ ppStack.config Index: branches/czw_branch/20170908/ippconfig/hsc/camera.config =================================================================== --- branches/czw_branch/20170908/ippconfig/hsc/camera.config (revision 40129) +++ branches/czw_branch/20170908/ippconfig/hsc/camera.config (revision 40179) @@ -140,4 +140,5 @@ z STR HSC-z y STR HSC-Y + NB816 STR NB0816 NB921 STR NB0921 END @@ -153,5 +154,5 @@ END -DVO.CAMERADIR STR megacam # Camera directory for DVO +DVO.CAMERADIR STR hsc # Camera directory for DVO # convert supplied FPA.OBSTYPE values to abstract exptype names @@ -170,5 +171,5 @@ RECIPES METADATA # Other recipes -# PSPHOT STR hsc/psphot.config # psphot details + PSPHOT STR hsc/psphot.config # psphot details PSASTRO STR hsc/psastro.config # psastro details # PSWARP STR hsc/pswarp.config # pswarp details Index: branches/czw_branch/20170908/ippconfig/hsc/format_cmf.config =================================================================== --- branches/czw_branch/20170908/ippconfig/hsc/format_cmf.config (revision 40129) +++ branches/czw_branch/20170908/ippconfig/hsc/format_cmf.config (revision 40179) @@ -194,4 +194,5 @@ # FPA.POSANGLE STR INR-STR FPA.POSANGLE STR INST-PA + FPA.ROTANGLE STR INST-PA FPA.RA STR RA FPA.DEC STR DEC Index: branches/czw_branch/20170908/ippconfig/hsc/format_mef.config =================================================================== --- branches/czw_branch/20170908/ippconfig/hsc/format_mef.config (revision 40129) +++ branches/czw_branch/20170908/ippconfig/hsc/format_mef.config (revision 40179) @@ -195,4 +195,5 @@ # FPA.POSANGLE STR INR-STR -- this is equivalent to the GPC1 ROT, not POS FPA.POSANGLE STR INST-PA + FPA.ROTANGLE STR INST-PA FPA.RA STR RA FPA.DEC STR DEC Index: branches/czw_branch/20170908/ippconfig/hsc/format_raw.config =================================================================== --- branches/czw_branch/20170908/ippconfig/hsc/format_raw.config (revision 40129) +++ branches/czw_branch/20170908/ippconfig/hsc/format_raw.config (revision 40179) @@ -194,4 +194,5 @@ # FPA.POSANGLE STR INR-STR FPA.POSANGLE STR INST-PA + FPA.ROTANGLE STR INST-PA FPA.RA STR RA FPA.DEC STR DEC Index: branches/czw_branch/20170908/ippconfig/hsc/ghost.model.mdc =================================================================== --- branches/czw_branch/20170908/ippconfig/hsc/ghost.model.mdc (revision 40179) +++ branches/czw_branch/20170908/ippconfig/hsc/ghost.model.mdc (revision 40179) @@ -0,0 +1,78 @@ +# The R_TERMS polynomials have been converted from um units to mm units to be consistent with C_TERMS polynomials. + +GHOST.MODEL.HSC MULTI + GHOST.MODEL.HSC METADATA + FILTER STR g + @C_TERMS F32 3.86129097e-04 -4.56480396e+02 5.94514717e+02 -8.20984903e+03 4.35455723e+04 -1.19834803e+05 1.60789252e+05 -8.66581086e+04 # a comment + @R_TERMS F32 2.80257250e+00 -1.14163334e+01 6.25378363e+02 -6.59305008e+03 3.21678893e+04 -7.68658356e+04 8.71449290e+04 -3.66203193e+04 # a comment + LIMIT F32 12.0 + SCALE F32 3.49794238e-06 +# WRONG ORDER +# @C_TERMS F32 -8.66581086e+04 1.60789252e+05 -1.19834803e+05 4.35455723e+04 -8.20984903e+03 5.94514717e+02 -4.56480396e+02 3.86129097e-04 +# @R_TERMS F32 -3.66203193e+04 8.71449290e+04 -7.68658356e+04 3.21678893e+04 -6.59305008e+03 6.25378363e+02 -1.14163334e+01 2.80257250e+00 + END + GHOST.MODEL.HSC METADATA + FILTER STR r + @C_TERMS F32 1.10217740e-03 -4.13994264e+02 -1.08162329e+02 -1.75333116e+02 -1.44358813e+03 7.03878935e+03 -1.55470840e+04 1.05167578e+04 # a comment + @R_TERMS F32 7.84632263e+00 -1.34007061e-02 9.57272709e+01 -9.45421731e+02 5.64254906e+03 -1.97396910e+04 3.55382290e+04 -2.46004825e+04 # a comment + LIMIT F32 11.5 + SCALE F32 3.49794238e-06 + # @C_TERMS for original r 8.14947896e+03 -1.24163720e+04 5.41143508e+03 -1.02024286e+03 -2.39667209e+02 -1.00805716e+02 -4.16757016e+02 -1.48799186e-03 + # @R_TERMS for original r -2.34342573e+07 3.47643396e+07 -1.98155551e+07 5.75428423e+06 -9.54233470e+05 9.14792638e+04 1.45742667e+03 6.96530721e+03 +# WRONG ORDER +# @C_TERMS F32 1.05167578e+04 -1.55470840e+04 7.03878935e+03 -1.44358813e+03 -1.75333116e+02 -1.08162329e+02 -4.13994264e+02 1.10217740e-03 +# @R_TERMS F32 -2.46004825e+04 3.55382290e+04 -1.97396910e+04 5.64254906e+03 -9.45421731e+02 9.57272709e+01 -1.34007061e-02 7.84632263e+00 + END + GHOST.MODEL.HSC METADATA + FILTER STR i + @C_TERMS F32 -3.95203567e-03 -4.55641159e+02 5.73358984e+02 -8.75512267e+03 5.43205903e+04 -1.74052246e+05 2.70933733e+05 -1.66554502e+05 # a comment + @R_TERMS F32 2.24351724e+00 -1.59383063e+01 8.56547254e+02 -1.09461298e+04 6.86646484e+04 -2.19342516e+05 3.46712176e+05 -2.15706075e+05 # a comment + LIMIT F32 11.0 + SCALE F32 3.49794238e-06 + # @C_TERMS for original i 7.09135701e+03 -1.08899009e+04 4.53311864e+03 -7.45070274e+02 -2.86196488e+02 -9.58537252e+01 -4.15997892e+02 -2.10038501e-03 + # @R_TERMS for original i -2.65938293e+07 3.98815563e+07 -2.30650423e+07 6.77205238e+06 -1.11878753e+06 1.03127466e+05 1.44010771e+03 6.83384243e+03 +# WRONG ORDER +# @C_TERMS F32 -1.66554502e+05 2.70933733e+05 -1.74052246e+05 5.43205903e+04 -8.75512267e+03 5.73358984e+02 -4.55641159e+02 -3.95203567e-03 +# @R_TERMS F32 -2.15706075e+05 3.46712176e+05 -2.19342516e+05 6.86646484e+04 -1.09461298e+04 8.56547254e+02 -1.59383063e+01 2.24351724e+00 + END + GHOST.MODEL.HSC METADATA + FILTER STR z + @C_TERMS F32 -1.79546522e-03 -4.13201008e+02 -1.08898155e+02 -1.51535018e+02 -1.51139264e+03 6.94362981e+03 -1.47396976e+04 9.59579373e+03 # a comment + @R_TERMS F32 7.39449641e+00 7.28959803e-01 1.06569049e+02 -1.18036333e+03 7.42160510e+03 -2.62797780e+04 4.70124953e+04 -3.23514338e+04 # a comment + LIMIT F32 11.0 + SCALE F32 3.49794238e-06 +# WRONG ORDER +# @C_TERMS F32 9.59579373e+03 -1.47396976e+04 6.94362981e+03 -1.51139264e+03 -1.51535018e+02 -1.08898155e+02 -4.13201008e+02 -1.79546522e-03 +# @R_TERMS F32 -3.23514338e+04 4.70124953e+04 -2.62797780e+04 7.42160510e+03 -1.18036333e+03 1.06569049e+02 7.28959803e-01 7.39449641e+00 + END + GHOST.MODEL.HSC METADATA + FILTER STR y + @C_TERMS F32 -2.04774887e-03 -4.50571358e+02 3.73884700e+02 -6.38727141e+03 4.28067607e+04 -1.47306815e+05 2.42465240e+05 -1.56137979e+05 # a comment + @R_TERMS F32 2.04967838e+00 -8.92675213e+00 6.03954493e+02 -7.97354793e+03 5.32635002e+04 -1.79071587e+05 2.94475198e+05 -1.88924528e+05 # a comment + LIMIT F32 11.0 + SCALE F32 3.49794238e-06 +# WRONG ORDER +# @C_TERMS F32 -1.56137979e+05 2.42465240e+05 -1.47306815e+05 4.28067607e+04 -6.38727141e+03 3.73884700e+02 -4.50571358e+02 -2.04774887e-03 +# @R_TERMS F32 -1.88924528e+05 2.94475198e+05 -1.79071587e+05 5.32635002e+04 -7.97354793e+03 6.03954493e+02 -8.92675213e+00 2.04967838e+00 + END + GHOST.MODEL.HSC METADATA + FILTER STR NB816 + @C_TERMS F32 -3.52735759e-03 -4.54174618e+02 5.11013724e+02 -8.05597593e+03 5.11417991e+04 -1.67422366e+05 2.65225596e+05 -1.65501517e+05 # a comment + @R_TERMS F32 2.17954931e+00 -1.39611991e+01 7.88077371e+02 -1.01750660e+04 6.49121728e+04 -2.10323885e+05 3.36261570e+05 -2.11111649e+05 # a comment + LIMIT F32 11.0 + SCALE F32 3.49794238e-06 +# WRONG ORDER +# @C_TERMS F32 -1.65501517e+05 2.65225596e+05 -1.67422366e+05 5.11417991e+04 -8.05597593e+03 5.11013724e+02 -4.54174618e+02 -3.52735759e-03 +# @R_TERMS F32 -2.11111649e+05 3.36261570e+05 -2.10323885e+05 6.49121728e+04 -1.01750660e+04 7.88077371e+02 -1.39611991e+01 2.17954931e+00 + END + GHOST.MODEL.HSC METADATA + FILTER STR NB921 + @C_TERMS F32 -2.99761113e-03 -4.52301392e+02 4.49346634e+02 -7.38661211e+03 4.83455130e+04 -1.62550708e+05 2.63017929e+05 -1.66995347e+05 # a comment + @R_TERMS F32 2.11290612e+00 -1.15228777e+01 7.02729746e+02 -9.20923423e+03 6.01493518e+04 -1.98625679e+05 3.22255104e+05 -2.04640222e+05 # a comment + LIMIT F32 11.0 + SCALE F32 3.49794238e-06 +# WRONG ORDER +# @C_TERMS F32 -1.66995347e+05 2.63017929e+05 -1.62550708e+05 4.83455130e+04 -7.38661211e+03 4.49346634e+02 -4.52301392e+02 -2.99761113e-03 +# @R_TERMS F32 -2.04640222e+05 3.22255104e+05 -1.98625679e+05 6.01493518e+04 -9.20923423e+03 7.02729746e+02 -1.15228777e+01 2.11290612e+00 + END +END Index: branches/czw_branch/20170908/ippconfig/hsc/ppStack.config =================================================================== --- branches/czw_branch/20170908/ippconfig/hsc/ppStack.config (revision 40129) +++ branches/czw_branch/20170908/ippconfig/hsc/ppStack.config (revision 40179) @@ -11,4 +11,5 @@ HSC-Y F32 0.0 NB0921 F32 0.0 + NB0816 F32 0.0 END @@ -22,4 +23,5 @@ HSC-Y F32 30.0 NB0921 F32 30.0 + NB0816 F32 30.0 END Index: branches/czw_branch/20170908/ippconfig/hsc/psastro.config =================================================================== --- branches/czw_branch/20170908/ippconfig/hsc/psastro.config (revision 40129) +++ branches/czw_branch/20170908/ippconfig/hsc/psastro.config (revision 40179) @@ -7,6 +7,6 @@ DVO.GETSTAR.PHOTCODE STR i DVO.GETSTAR.MAX.RHO F32 30000.0 -DVO.GETSTAR.MIN.MAG F32 12.0 -DVO.GETSTAR.MIN.MAG.INST F32 -15.0 +DVO.GETSTAR.MIN.MAG F32 -3.0 +DVO.GETSTAR.MIN.MAG.INST F32 -27.0 PHOTCODE.DATA MULTI @@ -39,4 +39,10 @@ ZEROPT F32 25.1 PHOTCODE STR y + GHOST_MAX_MAG F32 -25.0 +END +PHOTCODE.DATA METADATA + FILTER STR NB816 + ZEROPT F32 26.0 + PHOTCODE STR i GHOST_MAX_MAG F32 -25.0 END @@ -116,16 +122,119 @@ # Masking parameters -# REFSTAR_MASK BOOL TRUE -# REFSTAR_MASK_MAX_MAG F32 -15.0 -# REFSTAR_MASK_SATSTAR_MAG_MAX F32 -15.0 -# REFSTAR_MASK_SATSTAR_MAG_SLOPE F32 10.15 -# REFSTAR_MASK_SATSTAR_POS_ZERO F32 0.0 - -# REFSTAR_MASK_SATSPIKE_MAG_SLOPE F32 0.5188 -# REFSTAR_MASK_SATSPIKE_MAG_MAX F32 -14.644 -# REFSTAR_MASK_SATSPIKE_OFFSET F32 0.0 -# REFSTAR_MASK_SATSPIKE_WIDTH F32 6.0 -# REFSTAR_MASK_SATSPIKE_L0 F32 0 -# REFSTAR_MASK_SATSPIKE_WIDTH_SLOPE F32 0.0106 +REFSTAR_MASK BOOL TRUE +REFSTAR_MASK_MAX_MAG F32 -15.0 +REFSTAR_MASK_SATSTAR_MAG_MAX F32 -15.0 +REFSTAR_MASK_SATSTAR_MAG_SLOPE F32 4.45 # Was 10.15 +REFSTAR_MASK_SATSTAR_POS_ZERO F32 -15.0 + +REFSTAR_MASK_SATSPIKE_MAG_SLOPE F32 0.35 +REFSTAR_MASK_SATSPIKE_MAG_MAX F32 -14.5 +REFSTAR_MASK_SATSPIKE_OFFSET F32 0.0 +REFSTAR_MASK_SATSPIKE_WIDTH F32 6.0 +REFSTAR_MASK_SATSPIKE_L0 F32 0 +REFSTAR_MASK_SATSPIKE_WIDTH_SLOPE F32 0.0106 + +REFSTAR_MASK_BLEED BOOL TRUE + +REFSTAR_MASK_BLEED_MAG_MAX F32 -14.5 +REFSTAR_MASK_BLEED_MAG_SLOPE F32 2.0 + +REFSTAR_MASK_BLEED_ORIENTATION_X BOOL FALSE +REFSTAR_MASK_BLEED_MAG_MAX_X F32 -15.0 +REFSTAR_MASK_BLEED_MAG_SLOPE_X F32 5.0 + +REFSTAR_MASK_BLEED_ORIENTATION_Y BOOL TRUE +REFSTAR_MASK_BLEED_MAG_MAX_Y F32 -14.5 +REFSTAR_MASK_BLEED_MAG_SLOPE_Y F32 0.425 + +REFSTAR_MASK_GLINTS BOOL TRUE +GLINT_MAX_MAG F32 -19 +GLINT_LENGTH_MAG_SLOPE F32 0.0 +GLINT_LENGTH_MAG_ZERO F32 0.0 +GLINT_WIDTH F32 80.0 + + +# The slopes and radii and region below are all in FPA units, 1/13.5 of a regular pixel. +GLINT.REGION MULTI +GLINT.REGION METADATA + GLINT.TYPE STR HSC + GLINT.PIXEL.SCALE F32 3.49794238e-06 + REGION STR [-350000:350000,-350000:350000] + GLINT.RADIUS.INNER F32 252435 # 0.883 degrees + GLINT.RADIUS.OUTER F32 283025 # 0.990 degrees + + GLINT.FILTER.TERM MULTI + GLINT.FILTER.TERM METADATA + FILTER STR g + @C_TERMS F32 -4467.9759580157 10022.1277915657 -5156.3556611527 # a comment + @R_TERMS F32 -4375.9715493905 9717.6993803237 -5203.2424778605 # a comment + FACTOR F32 1.002 + LIMIT F32 14.0 + END + GLINT.FILTER.TERM METADATA + FILTER STR r + @C_TERMS F32 -4354.2575893763 9762.9766653532 -5011.1243975288 # a comment + @R_TERMS F32 -4263.2230039887 9457.4516871365 -5057.2762518308 # a comment + FACTOR F32 1.0015 + LIMIT F32 13.5 + END + GLINT.FILTER.TERM METADATA + FILTER STR i + @C_TERMS F32 -4243.1837126908 9517.3512135240 -4876.3237965836 # a comment + @R_TERMS F32 -4154.5095389959 9215.4140278639 -4924.2923196315 # a comment + FACTOR F32 1.001 + LIMIT F32 13.5 + END + GLINT.FILTER.TERM METADATA + FILTER STR z + @C_TERMS F32 -3966.2315539829 8920.7199277389 -4555.9567436784 # a comment + @R_TERMS F32 -3878.4529585585 8619.7241050446 -4604.3603386145 # a comment + FACTOR F32 0.9995 + LIMIT F32 13.0 + END + GLINT.FILTER.TERM METADATA + FILTER STR y + @C_TERMS F32 -3723.5919271089 8403.7641770854 -4281.0159171765 # a comment + @R_TERMS F32 -3637.0352789862 8104.8245491120 -4330.4583627761 # a comment + FACTOR F32 1.000 + LIMIT F32 13.0 + END + GLINT.FILTER.TERM METADATA + FILTER STR NB816 + @C_TERMS F32 -4243.1837126908 9517.3512135240 -4876.3237965836 # a comment + @R_TERMS F32 -4154.5095389959 9215.4140278639 -4924.2923196315 # a comment + FACTOR F32 1.001 + LIMIT F32 13.5 + END + GLINT.FILTER.TERM METADATA + FILTER STR NB921 + @C_TERMS F32 -3966.2315539829 8920.7199277389 -4555.9567436784 # a comment + @R_TERMS F32 -3878.4529585585 8619.7241050446 -4604.3603386145 # a comment + FACTOR F32 1.0005 + LIMIT F32 12.0 + END + +END + +# GLINT.RADIUS.INNER F32 245000 +# GLINT.RADIUS.OUTER F32 260000 + +# GLINT.ANGLE.INNER_SLOPE F32 -0.0059259 +# GLINT.ANGLE.INNER_ZERO F32 72 + +# GLINT.ANGLE.OUTER_SLOPE F32 0.0088889 +# GLINT.ANGLE.OUTER_ZERO F32 -41 + +# GLINT.ANGLE.INNER_SLOPE F32 0.0 +# GLINT.ANGLE.INNER_ZERO F32 7.0 +# GLINT.ANGLE.INNER_ZERO F32 1.5 +# GLINT.ANGLE.OUTER_SLOPE F32 0.0 +# GLINT.ANGLE.OUTER_ZERO F32 22.0 +# GLINT.ANGLE.OUTER_ZERO F32 6.0 + + + + + # REFSTAR_MASK_GLINTS BOOL TRUE @@ -137,2 +246,24 @@ # GLINT.TYPE STR HSC # END + +STATICSKY_CAL METADATA + REFSTAR_MASK BOOL FALSE + PSASTRO.SAVE.REFMATCH BOOL FALSE + PSASTRO.FIX.CHIPS BOOL FALSE + PSASTRO.USE.MODEL BOOL FALSE + PSASTRO.MOSAIC.MODE BOOL FALSE + DVO.GETSTAR.FIXED.ZEROPT F32 30.0 + DVO.GETSTAR.MIN.MAG.INST F32 -22.0 + PSASTRO.CATDIR STR PS1.REF.20140713 + ZERO.POINT.USE.MEAN BOOL TRUE + PSASTRO.MATCH.RADIUS.N0 F32 30 + PSASTRO.MATCH.RADIUS.N1 F32 20 + PSASTRO.MATCH.RADIUS.N2 F32 10 + PSASTRO.MATCH.RADIUS.N3 F32 2 + PSASTRO.MATCH.RADIUS.N4 F32 1 + PSASTRO.MATCH.FIT.NITER S32 5 + + PSASTRO.MAX.ERROR F32 0.2 # in arcsec units, if larger than 0.2" on skycell then a problem.. + + TRANSFORM.GALAXY.SHAPE.ERRORS BOOL T +END Index: branches/czw_branch/20170908/ippconfig/hsc/psphot.config =================================================================== --- branches/czw_branch/20170908/ippconfig/hsc/psphot.config (revision 40179) +++ branches/czw_branch/20170908/ippconfig/hsc/psphot.config (revision 40179) @@ -0,0 +1,88 @@ +BACKGROUND.XBIN S32 600 # size of background superpixels +BACKGROUND.YBIN S32 600 # size of background superpixels +BACKGROUND.XSAMPLE F32 2.0 # size of background superpixels +BACKGROUND.YSAMPLE F32 2.0 # size of background superpixels + +PSF_MODEL STR PS_MODEL_HSC_V1 + + + +STACKPHOT METADATA + SAVE.RESID BOOL TRUE + + PSF_MODEL STR PS_MODEL_HSC_V1 + EXTENDED_SOURCE_FITS BOOL FALSE # perform any of the aperture-like measurements? + EXTENDED_SOURCE_GALAXY_PARAMS BOOL FALSE # measure additional paramters for "galaxies" + EXTENDED_SOURCE_SN_LIM F32 10.0 + EXTENDED_SOURCE_PETROSIAN BOOL FALSE + EXTENDED_SOURCE_ANNULI BOOL TRUE + EXT.NSIGMA.LIMIT.USE BOOL FALSE + EXT.FIT.MIN.GAL.LIMIT.USE BOOL TRUE # limit extended source fits by galactic coordinates + + PSPHOT.STACK.MATCH.PSF.SOURCE STR AUTO # which inputs to convolve? (RAW, CNV, AUTO) + PSPHOT.STACK.TARGET.PSF.AUTO BOOL F # automatically determine target PSF size? + # NOTE: if PSPHOT.STACK.USE.RAW is set to False, KRON.SMOOTH should be + # set to false as well in order to get the best MOMENTS_R1 measurement + PSPHOT.STACK.USE.RAW BOOL T # perform photometry & morphology analysis on the raw image + KRON_SMOOTH BOOL T # smooth before kron radius measurement + @PSPHOT.STACK.TARGET.PSF.FWHM F32 6.0 8.0 # FWHM of target PSF (if NOT AUTO sized; pixels) + RADIAL_APERTURES BOOL T # calculate flux in circular radial apertures? + RADIAL_APERTURES_SN_LIM F32 0.0 # S/N limit for radial aperture calculation + OUTPUT.FORMAT STR PS1_SV4 + + # subtract radial profiles for extremely bright stars? + SUBTRACT_SATSTAR_PROFILE BOOL T + + CONSTANT_PHOTOMETRIC_WEIGHTS BOOL TRUE + POISSON.ERRORS.PHOT.LMM BOOL TRUE + + LMM_FIT_CHISQ_CONVERGENCE BOOL FALSE # use the (old) chisq convergence or the new parameter one + LMM_FIT_GAIN_FACTOR_MODE S32 2 # there are 3 options to determine the gain factor (0,1,2) + EXT_FIT_MIN_TOL F32 0.1 # Fit tolerance for EXT + + EXT_FIT_MAX_RADIUS F32 50.0 + PEAKS_NMAX_TOTAL S32 0 # set this to limit the allowed number of peaks - Yields fault instead of avoid memory explosion + + REPAIR.PETROSIAN.MAG.ERRORS BOOL FALSE # write PETRO_MAG_ERR as dFlux / Flux when writing + + + # PSPHOT.STACK.MATCH.FILTER + # list of metadata to control how source matching is performed for each band + # FILTER.ID is used to match to the FPA's concept value FPA.FILTERID + # ORDER controls the order in which the bands are processed. Lower numbers are processed first. + # The intent is to process in order of sensitivity + # MATCH.ALL if true causes forced photometry sources for objects detected in that band but don't meet + # the minimum number of detections cut + # TYPE is used in the recipe bookkeeping and has no semantic (but must be unique) + # Note: if an input has a filter not in this list it will be processed after all of the others with MATCH.ALL = FALSE + PSPHOT.STACK.MATCH.FILTERS METADATA + TYPE PSPHOT_STACK.MATCH.FILTER FILTER.ID ORDER MATCH.ALL Y.RATIO.MAX + gband PSPHOT_STACK.MATCH.FILTER HSC-g 4 TRUE 0 + rband PSPHOT_STACK.MATCH.FILTER HSC-r2 2 TRUE 0 + iband PSPHOT_STACK.MATCH.FILTER HSC-i2 1 TRUE 0 + zband PSPHOT_STACK.MATCH.FILTER HSC-z 3 TRUE 0 + yband PSPHOT_STACK.MATCH.FILTER HSC-Y 5 TRUE 0 + NB9 PSPHOT_STACK.MATCH.FILTER NB0921 6 TRUE 0 + NB8 PSPHOT_STACK.MATCH.FILTER NB0816 7 TRUE 0 + END + + EXT.ANALYSIS.MAG.LIMITS METADATA + TYPE DATA FILTER.ID MAG.LIMIT.PETRO MAG.LIMIT.EXTFIT + gband DATA g 28 25.5 + rband DATA r 28 25.5 + iband DATA i 27 24.5 + zband DATA z 26.5 23 + yband DATA y 26 22.5 + NB9 DATA NB0921 26 22.5 + NB8 DATA NB0816 26 22.5 + other DATA any 25 21.5 + +# gband DATA g 25 21.5 +# rband DATA r 25 21.5 +# iband DATA i 25 21.5 +# zband DATA z 25 20.5 +# yband DATA y 25 19.5 +# NB DATA NB0921 25 20.5 +# other DATA any 25 21.5 + END +END Index: branches/czw_branch/20170908/ippconfig/recipes/psastro.config =================================================================== --- branches/czw_branch/20170908/ippconfig/recipes/psastro.config (revision 40129) +++ branches/czw_branch/20170908/ippconfig/recipes/psastro.config (revision 40179) @@ -154,6 +154,16 @@ REFSTAR_MASK_BLEED BOOL FALSE + REFSTAR_MASK_BLEED_MAG_MAX F32 -15.0 REFSTAR_MASK_BLEED_MAG_SLOPE F32 5.0 + +# Additional parameters to allow for saturation bleeds to be masked in either direction. +REFSTAR_MASK_BLEED_ORIENTATION_X BOOL FALSE +REFSTAR_MASK_BLEED_MAG_MAX_X F32 -15.0 +REFSTAR_MASK_BLEED_MAG_SLOPE_X F32 5.0 + +REFSTAR_MASK_BLEED_ORIENTATION_Y BOOL FALSE +REFSTAR_MASK_BLEED_MAG_MAX_Y F32 -15.0 +REFSTAR_MASK_BLEED_MAG_SLOPE_Y F32 5.0 REFSTAR_MASK_GLINTS BOOL FALSE Index: branches/czw_branch/20170908/psModules/src/objects/Makefile.am =================================================================== --- branches/czw_branch/20170908/psModules/src/objects/Makefile.am (revision 40129) +++ branches/czw_branch/20170908/psModules/src/objects/Makefile.am (revision 40179) @@ -92,5 +92,6 @@ models/pmModel_EXP.c \ models/pmModel_DEV.c \ - models/pmModel_TRAIL.c + models/pmModel_TRAIL.c \ + models/pmModel_HSC_V1.c pkginclude_HEADERS = \ @@ -142,5 +143,6 @@ models/pmModel_EXP.h \ models/pmModel_DEV.h \ - models/pmModel_TRAIL.h + models/pmModel_TRAIL.h \ + models/pmModel_HSC_V1.h CLEANFILES = *~ Index: branches/czw_branch/20170908/psModules/src/objects/models/fwhm.sh =================================================================== --- branches/czw_branch/20170908/psModules/src/objects/models/fwhm.sh (revision 40129) +++ branches/czw_branch/20170908/psModules/src/objects/models/fwhm.sh (revision 40179) @@ -161,8 +161,8 @@ 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,15 +180,16 @@ 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,4 +207,66 @@ $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,5 +286,5 @@ 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,5 +292,5 @@ delete fwhm_v k_v - $SAVE = 0 + $SAVE = 1 if ($SAVE) exec rm -f $3 if ($SAVE) output $3 @@ -240,5 +303,5 @@ 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,8 +310,8 @@ 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 + Index: branches/czw_branch/20170908/psModules/src/objects/models/hscv1.fwhm.h =================================================================== --- branches/czw_branch/20170908/psModules/src/objects/models/hscv1.fwhm.h (revision 40179) +++ branches/czw_branch/20170908/psModules/src/objects/models/hscv1.fwhm.h (revision 40179) @@ -0,0 +1,111 @@ +# define FWHM_BIN 0.2 +# define MIN_FWHM_BIN -1 +# define N_FWHM_BIN 106 +static float HSC_V1_Scale[] = { +3.75, // -1.0, 1.76 +3.55, // -0.8, 1.57 +3.35, // -0.6, 1.40 +3.17, // -0.4, 1.25 +2.99, // -0.2, 1.12 +2.83, // 0.0, 1.00 +2.68, // 0.2, 0.90 +2.54, // 0.4, 0.81 +2.41, // 0.6, 0.73 +2.30, // 0.8, 0.66 +2.19, // 1.0, 0.60 +2.10, // 1.2, 0.55 +2.01, // 1.4, 0.51 +1.93, // 1.6, 0.47 +1.86, // 1.8, 0.43 +1.79, // 2.0, 0.40 +1.74, // 2.2, 0.38 +1.68, // 2.4, 0.35 +1.63, // 2.6, 0.33 +1.58, // 2.8, 0.31 +1.54, // 3.0, 0.30 +1.50, // 3.2, 0.28 +1.46, // 3.4, 0.27 +1.43, // 3.6, 0.25 +1.39, // 3.8, 0.24 +1.36, // 4.0, 0.23 +1.33, // 4.2, 0.22 +1.31, // 4.4, 0.21 +1.28, // 4.6, 0.20 +1.26, // 4.8, 0.20 +1.23, // 5.0, 0.19 +1.21, // 5.2, 0.18 +1.19, // 5.4, 0.18 +1.17, // 5.6, 0.17 +1.15, // 5.8, 0.17 +1.13, // 6.0, 0.16 +1.12, // 6.2, 0.16 +1.10, // 6.4, 0.15 +1.08, // 6.6, 0.15 +1.07, // 6.8, 0.14 +1.05, // 7.0, 0.14 +1.04, // 7.2, 0.14 +1.03, // 7.4, 0.13 +1.01, // 7.6, 0.13 +1.00, // 7.8, 0.13 +0.99, // 8.0, 0.12 +0.98, // 8.2, 0.12 +0.97, // 8.4, 0.12 +0.95, // 8.6, 0.11 +0.94, // 8.8, 0.11 +0.93, // 9.0, 0.11 +0.92, // 9.2, 0.11 +0.91, // 9.4, 0.10 +0.91, // 9.6, 0.10 +0.90, // 9.8, 0.10 +0.89, // 10.0, 0.10 +0.88, // 10.2, 0.10 +0.87, // 10.4, 0.09 +0.86, // 10.6, 0.09 +0.85, // 10.8, 0.09 +0.85, // 11.0, 0.09 +0.84, // 11.2, 0.09 +0.83, // 11.4, 0.09 +0.83, // 11.6, 0.09 +0.82, // 11.8, 0.08 +0.81, // 12.0, 0.08 +0.81, // 12.2, 0.08 +0.80, // 12.4, 0.08 +0.79, // 12.6, 0.08 +0.79, // 12.8, 0.08 +0.78, // 13.0, 0.08 +0.77, // 13.2, 0.08 +0.77, // 13.4, 0.07 +0.76, // 13.6, 0.07 +0.76, // 13.8, 0.07 +0.75, // 14.0, 0.07 +0.75, // 14.2, 0.07 +0.74, // 14.4, 0.07 +0.74, // 14.6, 0.07 +0.73, // 14.8, 0.07 +0.73, // 15.0, 0.07 +0.72, // 15.2, 0.07 +0.72, // 15.4, 0.06 +0.71, // 15.6, 0.06 +0.71, // 15.8, 0.06 +0.70, // 16.0, 0.06 +0.70, // 16.2, 0.06 +0.70, // 16.4, 0.06 +0.69, // 16.6, 0.06 +0.69, // 16.8, 0.06 +0.68, // 17.0, 0.06 +0.68, // 17.2, 0.06 +0.68, // 17.4, 0.06 +0.67, // 17.6, 0.06 +0.67, // 17.8, 0.06 +0.66, // 18.0, 0.06 +0.66, // 18.2, 0.05 +0.66, // 18.4, 0.05 +0.65, // 18.6, 0.05 +0.65, // 18.8, 0.05 +0.65, // 19.0, 0.05 +0.64, // 19.2, 0.05 +0.64, // 19.4, 0.05 +0.64, // 19.6, 0.05 +0.63, // 19.8, 0.05 +0.63, // 20.0, 0.05 +}; Index: branches/czw_branch/20170908/psModules/src/objects/models/pmModel_HSC_V1.c =================================================================== --- branches/czw_branch/20170908/psModules/src/objects/models/pmModel_HSC_V1.c (revision 40179) +++ branches/czw_branch/20170908/psModules/src/objects/models/pmModel_HSC_V1.c (revision 40179) @@ -0,0 +1,505 @@ +/****************************************************************************** + * this file defines the HSC_V1 source shape model. Note that these model functions are + * loaded by pmModelClass.c using 'include', and thus need no 'include' statements of + * their own. The models use a psVector to represent the set of parameters, with the + * sequence used to specify the meaning of the parameter. The meaning of the parameters + * may thus vary depending on the specifics of the model. All models which are used as a + * PSF representations share a few parameters, for which # define names are listed in + * pmModel.h: + + power-law with fitted linear term + 1 / (1 + kz + z^(3.33/2)) (z = r^2, so r^3.33) + + * PM_PAR_SKY 0 - local sky : note that this is unused and may be dropped in the future + * PM_PAR_I0 1 - central intensity + * PM_PAR_XPOS 2 - X center of object + * PM_PAR_YPOS 3 - Y center of object + * PM_PAR_SXX 4 - X^2 term of elliptical contour (SigmaX / sqrt(2)) + * PM_PAR_SYY 5 - Y^2 term of elliptical contour (SigmaY / sqrt(2)) + * PM_PAR_SXY 6 - X*Y term of elliptical contour + * PM_PAR_7 7 - amplitude of the linear component (k) + *****************************************************************************/ + +#include +#include +#include "pmHDU.h" +#include "pmFPA.h" + +#include "pmTrend2D.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmSpan.h" +#include "pmFootprintSpans.h" +#include "pmFootprint.h" +#include "pmPeaks.h" +#include "pmMoments.h" +#include "pmModelFuncs.h" +#include "pmModelClass.h" +#include "pmModel.h" +#include "pmModelUtils.h" +#include "pmSourceMasks.h" +#include "pmSourceExtendedPars.h" +#include "pmSourceDiffStats.h" +#include "pmSourceSatstar.h" +#include "pmSourceLensing.h" +#include "pmSource.h" +#include "pmSourceFitModel.h" +#include "pmPSF.h" +#include "pmPSFtry.h" +#include "pmDetections.h" + +#include "pmModel_HSC_V1.h" + +# define PM_MODEL_NPARAM 8 +# define PM_MODEL_FUNC pmModelFunc_HSC_V1 +# define PM_MODEL_FLUX pmModelFlux_HSC_V1 +# define PM_MODEL_GUESS pmModelGuess_HSC_V1 +# define PM_MODEL_LIMITS pmModelLimits_HSC_V1 +# define PM_MODEL_RADIUS pmModelRadius_HSC_V1 +# define PM_MODEL_SET_FWHM pmModelSetFWHM_HSC_V1 +# define PM_MODEL_FROM_PSF pmModelFromPSF_HSC_V1 +# define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_HSC_V1 +# define PM_MODEL_FIT_STATUS pmModelFitStatus_HSC_V1 +# define PM_MODEL_SET_LIMITS pmModelSetLimits_HSC_V1 + +# define ALPHA 1.8 +# define ALPHA_M 0.8 + +// the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) +// 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords +// values need to be pixel coords + +// Lax parameter limits +static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 }; +static float paramsMaxLax[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; + +// Moderate parameter limits +// Tolerate a small divot (k < 0) +static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 }; +static float paramsMaxModerate[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; + +// Strict parameter limits +// k = PAR_7 < 0 is very undesirable (big divot in the middle) +static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 }; +static float paramsMaxStrict[] = { 1.0e5, 1.0e9, 1.0e5, 1.0e5, 100, 100, 1.0, 20.0 }; + +// Parameter limits to use +static float *paramsMinUse = paramsMinLax; +static float *paramsMaxUse = paramsMaxLax; +static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 }; + +static bool limitsApply = true; // Apply limits? + +psF32 PM_MODEL_FUNC (psVector *deriv, + const psVector *params, + const psVector *pixcoord) +{ + psF32 *PAR = params->data.F32; + + psF32 X = pixcoord->data.F32[0] - PAR[PM_PAR_XPOS]; + psF32 Y = pixcoord->data.F32[1] - PAR[PM_PAR_YPOS]; + psF32 px = X / PAR[PM_PAR_SXX]; + psF32 py = Y / PAR[PM_PAR_SYY]; + psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; + + // XXX if the elliptical contour is defined in valid way, this step should not be required. + // other models (like PGAUSS) don't use fractional powers, and thus do not have NaN values + // for negative values of z + // XXX use an assert here to force the elliptical parameters to be correctly determined + // if (z < 0) z = 0; + assert (z >= 0); + + psF32 zp = pow(z,ALPHA_M); + psF32 r = 1.0 / (1 + PAR[PM_PAR_7]*z + z*zp); + + psF32 r1 = PAR[PM_PAR_I0]*r; + psF32 f = r1 + PAR[PM_PAR_SKY]; + + if (deriv != NULL) { + psF32 *dPAR = deriv->data.F32; + + // note difference from a pure gaussian: q = params->data.F32[PM_PAR_I0]*r + psF32 t = r1*r; + psF32 q = t*(PAR[PM_PAR_7] + ALPHA*zp); + + dPAR[PM_PAR_SKY] = +1.0; + dPAR[PM_PAR_I0] = +r; + dPAR[PM_PAR_XPOS] = q*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); + dPAR[PM_PAR_YPOS] = q*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); + // the extra factor of 2 below is needed to avoid excessive swings + dPAR[PM_PAR_SXX] = +4.0*q*px*px/PAR[PM_PAR_SXX]; + dPAR[PM_PAR_SYY] = +4.0*q*py*py/PAR[PM_PAR_SYY]; + dPAR[PM_PAR_SXY] = -q*X*Y; + dPAR[PM_PAR_7] = -t*z; + } + return(f); +} + +// define the parameter limits +// AR_MAX is the maximum allowed axis ratio +// AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2) +# define AR_MAX 20.0 +# define AR_RATIO 0.99 +bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) +{ + if (!limitsApply) { + return true; + } + psAssert(nParam >= 0 && nParam < PM_MODEL_NPARAM, "Parameter index is out of bounds"); + + // we need to calculate the limits for SXY specially + float q2 = NAN; + if (nParam == PM_PAR_SXY) { + float f1 = 2.0 / PS_SQR(params[PM_PAR_SYY]) + 2.0 / PS_SQR(params[PM_PAR_SXX]); + float f2 = 2.0 / PS_SQR(params[PM_PAR_SYY]) - 2.0 / PS_SQR(params[PM_PAR_SXX]); + float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); + q1 = (q1 < 0.0) ? 0.0 : q1; + // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg.. Saturate at that + // angle and let f2,f1 fight it out + q2 = 0.5*sqrtf(q1); + } + + switch (mode) { + case PS_MINIMIZE_BETA_LIMIT: { + psAssert(beta, "Require beta to limit beta"); + float limit = betaUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (fabs(beta[nParam]) > fabs(limit)) { + beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit); + psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", + nParam, beta[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MIN: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMinUse, "Require parameter limits to limit parameters"); + float limit = paramsMinUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] < limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + case PS_MINIMIZE_PARAM_MAX: { + psAssert(params, "Require parameters to limit parameters"); + psAssert(paramsMaxUse, "Require parameter limits to limit parameters"); + float limit = paramsMaxUse[nParam]; + if (nParam == PM_PAR_SXY) { + limit *= q2; + } + if (params[nParam] > limit) { + params[nParam] = limit; + psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", + nParam, params[nParam], limit); + return false; + } + return true; + } + default: + psAbort("invalid choice for limits"); + } + psAbort("should not reach here"); + return false; +} + + +// make an initial guess for parameters +// 0.5 PIX: moments and peaks are in pixel coords, thus so are model parameters +bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal) +{ + psF32 *PAR = model->params->data.F32; + + // sky is set to 0.0 + PAR[PM_PAR_SKY] = 0.0; + + // set the shape parameters + if (!pmModelSetShape(&PAR[PM_PAR_SXX], &PAR[PM_PAR_SXY], &PAR[PM_PAR_SYY], source->moments, false)) { + return false; + } + + // set the model normalization + if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { + return false; + } + + // set the model position + if (!pmModelSetPosition(&PAR[PM_PAR_XPOS], &PAR[PM_PAR_YPOS], source)) { + return false; + } + + // extra parameter + PAR[PM_PAR_7] = 0.5; + + return(true); +} + +// integrate the model to get the full flux +psF64 PM_MODEL_FLUX (const psVector *params) +{ + float z, norm; + + psF32 *PAR = params->data.F32; + + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); + float AspectRatio = axes.minor / axes.major; + + // flux = 2 \pi \int f(r) r dr + norm = 0.0; + + # define DR 0.25 + + // f = f(r) * r + float f0 = 0.0; + float f1, f2; + for (float r = DR; r < 150; r += DR) { + z = 0.5 * PS_SQR(r / axes.major); + f1 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); + r += DR; + z = 0.5 * PS_SQR(r / axes.major); + f2 = r / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); + norm += f0 + 4*f1 + f2; + f0 = f2; + } + norm *= DR / 3.0; + + psF64 Flux = PAR[PM_PAR_I0] * norm * 2.0 * M_PI * AspectRatio; + + return(Flux); +} + +// define this function so it never returns Inf or NaN +// return the radius which yields the requested flux +psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) +{ + psF64 z; + + psF32 *PAR = params->data.F32; + + if (flux <= 0) return 1.0; + if (PAR[PM_PAR_I0] <= 0) return 1.0; + if (flux >= PAR[PM_PAR_I0]) return 1.0; + + psEllipseAxes axes; + pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], false); + psF64 sigma = axes.major; + + // we can solve directly if PAR_7 goes to 0.0 + if (PAR[PM_PAR_7] == 0.0) { + psF32 z = powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA); + return ( sigma * sqrt (2.0 * z) ); + } + psF64 limit = flux / PAR[PM_PAR_I0]; + + // use the fact that f is monotonically decreasing + z = 0; + + // choose a z value guaranteed to be beyond our limit + float z0 = 0.0; + float z1 = pow((1.0 / limit), (1.0 / ALPHA)); + psAssert (isfinite(z1), "fix this code: z1 should not be nan for Io = %f, flux = %f, major = %f (%f, %f, %f)", + PAR[PM_PAR_I0], flux, axes.major, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY]); + if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0; + + // starting guess: + z = 0.5*(z0 + z1); + float dz = 1.0; + + // use Newton-Raphson to minimize f(z) - limit = 0 + for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) { + float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA)); + float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0)); + + float f = 1.0 / q; + float dfdz = -dqdz * f / q; + + dz = (f - limit) / dfdz; + + // fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q); + z -= dz; + z = PS_MAX(z, 0.0); + } + psF64 radius = sigma * sqrt (2.0 * z); + + if (isnan(radius)) + psAbort("error in code: radius is NaN"); + + return (radius); +} + +// I used the script in models/fwhm.sh to generate the trend of FWHM scaling vs the K value: +# include "hscv1.fwhm.h" + +psF64 PM_MODEL_SET_FWHM (const psVector *params, psF64 sigma) { + + psF32 *PAR = params->data.F32; + + float core = PAR[PM_PAR_7]; + + if (!isfinite(core)) return (2.0*M_SQRT2*sigma); + + int binCore = MAX(0, MIN (N_FWHM_BIN - 2, (int)((core - MIN_FWHM_BIN)/FWHM_BIN))); + float coreBin = binCore * FWHM_BIN + MIN_FWHM_BIN; + + float scale = (core - coreBin) * (HSC_V1_Scale[binCore + 1] - HSC_V1_Scale[binCore]) / FWHM_BIN + HSC_V1_Scale[binCore]; + + // fprintf (stderr, "core: %f, binCore: %d, scale: %f\n", core, binCore, scale); + + return (scale * sigma); +} + +bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) +{ + psF32 *out = modelPSF->params->data.F32; + psF32 *in = modelFLT->params->data.F32; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + for (int i = 0; i < psf->params->n; i++) { + if (psf->params->data[i] == NULL) { + out[i] = in[i]; + } else { + pmTrend2D *trend = psf->params->data[i]; + out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + } + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + bool useReff = modelPSF->class->useReff; + if (!pmPSF_FitToModel (out, 0.1, useReff)) { + psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]); + return false; + } + + // apply the model limits here: this truncates excessive extrapolation + // XXX do we need to do this still? should we put in asserts to test? + for (int i = 0; i < psf->params->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->data[i] == NULL) + continue; + + bool status = true; + status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MIN, i, out, NULL); + status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MAX, i, out, NULL); + if (!status) { + psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", + in[PM_PAR_XPOS], in[PM_PAR_YPOS]); + modelPSF->flags |= PM_MODEL_STATUS_LIMITS; + } + } + + return true; +} + +// construct the PSF model from the FLT model and the psf +// XXX is this sufficiently general do be a global function, not a pmModelClass function? +bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io) +{ + psF32 *PAR = model->params->data.F32; + + // we require these two parameters to exist + assert (psf->params->n > PM_PAR_YPOS); + assert (psf->params->n > PM_PAR_XPOS); + + PAR[PM_PAR_SKY] = 0.0; + PAR[PM_PAR_I0] = Io; + PAR[PM_PAR_XPOS] = Xo; + PAR[PM_PAR_YPOS] = Yo; + + // supply the model-fitted parameters, or copy from the input + for (int i = 0; i < psf->params->n; i++) { + if (i == PM_PAR_SKY) continue; + if (i == PM_PAR_I0) continue; + if (i == PM_PAR_XPOS) continue; + if (i == PM_PAR_YPOS) continue; + pmTrend2D *trend = psf->params->data[i]; + PAR[i] = pmTrend2DEval(trend, Xo, Yo); + } + + // the 2D PSF model fits polarization terms (E0,E1,E2) + // convert to shape terms (SXX,SYY,SXY) + // XXX user-defined value for limit? + bool useReff = model->class->useReff; + if (!pmPSF_FitToModel (PAR, 0.1, useReff)) { + psTrace ("psModules.objects", 3, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); + return false; + } + + // apply the model limits here: this truncates excessive extrapolation + // XXX do we need to do this still? should we put in asserts to test? + for (int i = 0; i < psf->params->n; i++) { + // apply the limits to all components or just the psf-model parameters? + if (psf->params->data[i] == NULL) + continue; + + bool status = true; + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL); + status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL); + if (!status) { + psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo); + model->flags |= PM_MODEL_STATUS_LIMITS; + } + } + return(true); +} + +bool PM_MODEL_FIT_STATUS (pmModel *model) +{ + bool status; + + psF32 *PAR = model->params->data.F32; + psF32 *dPAR = model->dparams->data.F32; + + status = true; + status &= (PAR[PM_PAR_I0] > 0); + status &= ((dPAR[PM_PAR_I0]/PAR[PM_PAR_I0]) < 0.5); + + return status; +} + + +void PM_MODEL_SET_LIMITS(pmModelLimitsType type) +{ + switch (type) { + case PM_MODEL_LIMITS_NONE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = true; + break; + case PM_MODEL_LIMITS_IGNORE: + paramsMinUse = NULL; + paramsMaxUse = NULL; + limitsApply = false; + break; + case PM_MODEL_LIMITS_LAX: + paramsMinUse = paramsMinLax; + paramsMaxUse = paramsMaxLax; + limitsApply = true; + break; + case PM_MODEL_LIMITS_MODERATE: + paramsMinUse = paramsMinModerate; + paramsMaxUse = paramsMaxModerate; + limitsApply = true; + break; + case PM_MODEL_LIMITS_STRICT: + paramsMinUse = paramsMinStrict; + paramsMaxUse = paramsMaxStrict; + limitsApply = true; + break; + default: + psAbort("Unrecognised model limits type: %x", type); + } + return; +} Index: branches/czw_branch/20170908/psModules/src/objects/models/pmModel_HSC_V1.h =================================================================== --- branches/czw_branch/20170908/psModules/src/objects/models/pmModel_HSC_V1.h (revision 40179) +++ branches/czw_branch/20170908/psModules/src/objects/models/pmModel_HSC_V1.h (revision 40179) @@ -0,0 +1,16 @@ +#ifndef PM_MODEL_HSC_V1_H + +#include "pmModel.h" + +psF32 pmModelFunc_HSC_V1(psVector *deriv, const psVector *params, const psVector *pixcoord); +bool pmModelLimits_HSC_V1(psMinConstraintMode mode, int nParam, float *params, float *beta); +bool pmModelGuess_HSC_V1(pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal); +psF64 pmModelFlux_HSC_V1(const psVector *params); +psF64 pmModelRadius_HSC_V1(const psVector *params, psF64 flux); +psF64 pmModelSetFWHM_HSC_V1(const psVector *params, psF64 flux); +bool pmModelFromPSF_HSC_V1(pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf); +bool pmModelParamsFromPSF_HSC_V1(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); +bool pmModelFitStatus_HSC_V1(pmModel *model); +void pmModelSetLimits_HSC_V1(pmModelLimitsType type); + +#endif Index: branches/czw_branch/20170908/psModules/src/objects/pmModelClass.c =================================================================== --- branches/czw_branch/20170908/psModules/src/objects/pmModelClass.c (revision 40129) +++ branches/czw_branch/20170908/psModules/src/objects/pmModelClass.c (revision 40179) @@ -47,4 +47,5 @@ # 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,4 +64,5 @@ {"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 }, }; Index: branches/czw_branch/20170908/psModules/src/objects/pmSourceOutputs.c =================================================================== --- branches/czw_branch/20170908/psModules/src/objects/pmSourceOutputs.c (revision 40129) +++ branches/czw_branch/20170908/psModules/src/objects/pmSourceOutputs.c (revision 40179) @@ -135,4 +135,7 @@ 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); Index: branches/czw_branch/20170908/psastro/src/psastroLoadGhosts.c =================================================================== --- branches/czw_branch/20170908/psastro/src/psastroLoadGhosts.c (revision 40129) +++ branches/czw_branch/20170908/psastro/src/psastroLoadGhosts.c (revision 40179) @@ -40,4 +40,7 @@ * calculate ghost FPA and Chip positions for the stars loaded on the FPA */ + + + bool psastroLoadGhosts (pmConfig *config) { @@ -256,2 +259,233 @@ return false; } + + + +bool psastroLoadGhostsHSC (pmConfig *config) { + bool status; + pmChip *chip = NULL; + pmCell *cell = NULL; + pmReadout *readout = NULL; + float zeropt, exptime, MAX_MAG; + psMetadata *ghostModel = NULL; + psVector *C_terms; + psVector *R_terms; + float glintPixelScale; + + psLogMsg ("psastro", PS_LOG_INFO, "determine ghost positions"); + + // select the current recipe + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSASTRO_RECIPE); + if (!recipe) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe"); + return false; + } + + bool REFSTAR_MASK_GHOST = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_GHOST"); + if (!REFSTAR_MASK_GHOST) return true; + + // Load model from file + char *ghostFile = psMetadataLookupStr (&status, recipe, "GHOST_MODEL"); + if (!strcasecmp(ghostFile, "NONE")) return true; + + if (!pmConfigFileRead (&ghostModel, ghostFile, "GHOST MODEL")) { + psError(PSASTRO_ERR_CONFIG, true, "Trouble loading ghost model"); + return false; + } + + // Allocate FPA and refstars + pmFPAview *view = pmFPAviewAlloc (0); + // select the input astrometry data (also carries the refstars) + pmFPAfile *astrom = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); + if (!astrom) { + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data"); + goto escape; + } + pmFPA *fpa = astrom->fpa; + + char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID"); + psMetadataItem *item = psMetadataLookup(glintItem->data.md, "GLINT.FILTER.TERM"); + if (!item) { + psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM data missing"); + return false; + } + if (item->type != PS_DATA_METADATA_MULTI) { + psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM not multi"); + return false; + } + psListIterator *iter = psListIteratorAlloc(item->data.list,PS_LIST_HEAD, false); + psMetadataItem *refItem = NULL; + while ((refItem = psListGetAndIncrement(iter))) { + if (refItem->type != PS_DATA_METADATA) { + psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM entry not metadata"); + return false; + } + char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); + if (!status) { + // psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + if (strcmp(refFilter, filter)) continue; + + C_terms = psMetadataLookupPtr(&status, refItem->data.md, "C_TERMS"); + R_terms = psMetadataLookupPtr(&status, refItem->data.md, "R_TERMS"); + MAX_MAG = psMetadataLookupF32(&status, refItem->data.md, "LIMIT"); + glintPixelScale = psMetadataLookupF32(&status, refItem->data.md, "SCALE"); + } + + // really error-out here? or just skip? + // if (!psastroZeroPointFromRecipe (&zeropt, &exptime, &MAX_MAG, fpa, recipe)) { + // psLogMsg ("psastro", PS_LOG_INFO, "failed to load zeropt data from recipe"); + // goto escape; + // } + + // recipe values are given in instrumental magnitudes + // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime) + // float MagOffset = zeropt + 2.5*log10(exptime); + // MAX_MAG += MagOffset; + + // this loop selects the matched stars for all chips + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (!chip->process || !chip->file_exists) { continue; } + if (!chip->fromFPA) { continue; } + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { continue; } + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } + + // select the raw objects for this readout (loaded in psastroChooseRefstars.c) + // XXX : note that we place limits on the refstar sample in psastroChooseRefstars.c: + // 1) on chip and 2) < PSASTRO.MAX.NREF. magnitude limits and clump exclusion are only + psArray *refstars = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.REFSTARS"); + if (refstars == NULL) { continue; } + + // identify the bright stars of interest + for (int i = 0; i < refstars->n; i++) { + pmAstromObj *ref = refstars->data[i]; + psTrace("psastro.ghost",5,"Begin ghost %d/%ld: MAX_MAG: %g; ref_mag: %g @ (%.10g,%.10g) TPA: (%f %f)", + i,refstars->n,MAX_MAG,ref->Mag,ref->sky->r * 180 / M_PI,ref->sky->d * 180.0 / M_PI,ref->TP->x,ref->RP->y); + if (ref->Mag > MAX_MAG) continue; + + psastroGhost *ghost = psastroGhostAlloc (); + + double R_TPA = sqrt(pow(star->TP->y,2) + pow(star->TP->x,2)); + double theta0 = atan2(star->TP->y,star->TP->x); + double star_radius_deg = R_TPA * glintPixelScale; + if (star_radius_deg > 0.5) { continue; } + + // Calculate expected position. + double C = C_terms->data.F32[0] + C_terms->data.F32[1] * star_radius_deg + C_terms->data.F32[2] * pow(star_radius_deg,2) + + C_terms->data.F32[3] * pow(star_radius_deg,3) + C_terms->data.F32[4] * pow(star_radius_deg,4) + + C_terms->data.F32[5] * pow(star_radius_deg,5) + C_terms->data.F32[6] * pow(star_radius_deg,6) + + C_terms->data.F32[7] * pow(star_radius_deg,7); + double R = R_terms->data.F32[0] + R_terms->data.F32[1] * star_radius_deg + R_terms->data.F32[2] * pow(star_radius_deg,2) + + R_terms->data.F32[3] * pow(star_radius_deg,3) + R_terms->data.F32[4] * pow(star_radius_deg,4) + + R_terms->data.F32[5] * pow(star_radius_deg,5) + R_terms->data.F32[6] * pow(star_radius_deg,6) + + R_terms->data.F32[7] * pow(star_radius_deg,7); + R = 13.5 * (R / 0.015); // Convert to TPA units. + + psPlane *fp = psPlaneAlloc(); + psPlane *tp = psPlaneAlloc(); + + tp->x = 13.5 * (C * cos(theta0) / 0.015 + 12.78); + tp->y = 13.5 * (C * sin(theta0) / 0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + + ghost->srcFP->x = ref->FP->x; + ghost->srcFP->y = ref->FP->y; + ghost->FP->x = fp->x; + ghost->FP->y = fp->y; + + ghost->inner.major = 0.0; + ghost->inner.minor = 0.0; + ghost->outer.major = R; + ghost->outer.minor = R; + + pmChip *ghostChip = psastroFindChip (&ghost->chip->x, &ghost->chip->y, fpa, ghost->FP->x, ghost->FP->y); + + + + if (ghostChip) { + psTrace("psastro.ghost",5," in ghost: %d/%ld: ref: (%f,%f) or (%s) ghost: (%f,%f) or (%f,%f,%s) %f inner: (%f,%f,%f) outer: (%f,%f,%f)", + i,refstars->n, + ref->FP->x,ref->FP->y, + psMetadataLookupStr(NULL,chip->concepts,"CHIP.NAME"), + ghost->FP->x,ghost->FP->y, + ghost->chip->x,ghost->chip->y,psMetadataLookupStr(NULL,ghostChip->concepts,"CHIP.NAME"), + rSrc, + ghost->inner.major,ghost->inner.minor,ghost->inner.theta, + ghost->outer.major,ghost->outer.minor,ghost->outer.theta); +/* psWarning(" GHOST_DATA: %d/%ld: ref: (%f,%f) or (%f,%f,%s) ghost: (%f,%f) or (%f,%f,%s) %f inner: (%f,%f,%f) outer: (%f,%f,%f)", */ +/* i,refstars->n, */ +/* ref->FP->x,ref->FP->y, */ +/* ref_chip_x,ref_chip_y,psMetadataLookupStr(NULL,chip->concepts,"CHIP.NAME"), */ +/* ghost->FP->x,ghost->FP->y, */ +/* ghost->chip->x,ghost->chip->y,psMetadataLookupStr(NULL,ghostChip->concepts,"CHIP.NAME"), */ +/* rSrc, */ +/* ghost->inner.major,ghost->inner.minor,ghost->inner.theta, */ +/* ghost->outer.major,ghost->outer.minor,ghost->outer.theta); */ + } + else { + psTrace("psastro.ghost",5," in ghost: %d/%ld: ref: (%f,%f) ghost: (%f,%f) or (%f,%f,%s) inner: (%f,%f,%f) outer: (%f,%f,%f)", + i,refstars->n, + ref->FP->x,ref->FP->y, + ghost->FP->x,ghost->FP->y, + ghost->chip->x,ghost->chip->y,"NONE", + ghost->inner.major,ghost->inner.minor,ghost->inner.theta, + ghost->outer.major,ghost->outer.minor,ghost->outer.theta); + } + + + + if (!ghostChip) goto skip; + if (!ghostChip->cells) goto skip; + if (!ghostChip->cells->n) goto skip; + + + pmCell *ghostCell = ghostChip->cells->data[0]; + if (!ghostCell) goto skip; + if (!ghostCell->readouts) goto skip; + if (!ghostCell->readouts->n) goto skip; + pmReadout *ghostReadout = ghostCell->readouts->data[0]; + if (!ghostReadout) goto skip; + + psArray *ghosts = psMetadataLookupPtr (&status, ghostReadout->analysis, "PSASTRO.GHOSTS"); + if (ghosts == NULL) { + ghosts = psArrayAllocEmpty (100); + if (!psMetadataAdd (ghostReadout->analysis, PS_LIST_TAIL, "PSASTRO.GHOSTS", PS_DATA_ARRAY, "astrometry matches", ghosts)) { + psError(PSASTRO_ERR_CONFIG, false, "failure to add ghosts to readout"); + goto escape; + } + psFree (ghosts); + } + + psArrayAdd (ghosts, 100, ghost); + + skip: + + psFree (ghost); + } + } + } + } + + psastroExtractFreeChipBounds(); + + psFree (C_terms); + psFree (R_terms); + psFree (ghostModel); + psFree (view); + return true; + +escape: + psFree (C_terms); + psFree (R_terms); + psFree (ghostModel); + psFree (view); + return false; +} Index: branches/czw_branch/20170908/psastro/src/psastroLoadGlints.c =================================================================== --- branches/czw_branch/20170908/psastro/src/psastroLoadGlints.c (revision 40129) +++ branches/czw_branch/20170908/psastro/src/psastroLoadGlints.c (revision 40179) @@ -44,4 +44,6 @@ double GLINT_LENGTH_MAG_ZERO = psMetadataLookupF32 (&status, recipe, "GLINT_LENGTH_MAG_ZERO"); double glintWidth = psMetadataLookupF32 (&status, recipe, "GLINT_WIDTH"); + double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE"); + // select the set of glint regions (GLINT.REGION is a MULTI of METADATA items) @@ -95,5 +97,5 @@ psProject (star->TP, star->sky, fpa->toSky); psPlaneTransformApply (star->FP, fpa->fromTPA, star->TP); - fprintf (stderr, "glint: %7.2f @ %8.1f, %8.1f\n", star->Mag, star->FP->x, star->FP->y); + fprintf (stderr, "glint: %7.2f @ %8.1f, %8.1f (%f %f) %8.1f %8.1f\n", star->Mag, star->FP->x, star->FP->y, star->sky->r * PS_DEG_RAD, star->sky->d * PS_DEG_RAD, star->TP->x,star->TP->y); // find the GLINT.REGION this star lands in (if any) @@ -298,4 +300,230 @@ } } + if (!strcasecmp(glintType, "HSC")) { + // It's inefficient to keep looking these up. + double GLINT_RADIUS_INNER = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.RADIUS.INNER"); + double GLINT_RADIUS_OUTER = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.RADIUS.OUTER"); + + // double R_FPA = sqrt(pow(star->FP->y,2) + pow(star->FP->x,2)); + double R_FPA = sqrt(pow(star->TP->y,2) + pow(star->TP->x,2)); + // R_FPA *= 1.30; + if (R_FPA < GLINT_RADIUS_INNER) { continue; } + if (R_FPA > GLINT_RADIUS_OUTER) { continue; } + psTrace("psastro.masks",4,"HSC_GLINT_STARS: %f %f %f\n", + // star->FP->x,star->FP->y,star->Mag); + star->TP->x,star->TP->y,star->Mag); + psVector *C_terms; + psVector *R_terms; + double scale_factor = 1.0; + + char *filter = psMetadataLookupStr (&status, fpa->concepts, "FPA.FILTERID"); + psMetadataItem *item = psMetadataLookup(glintItem->data.md, "GLINT.FILTER.TERM"); + if (!item) { + psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM data missing"); + return false; + } + if (item->type != PS_DATA_METADATA_MULTI) { + psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM not multi"); + return false; + } + psListIterator *iter = psListIteratorAlloc(item->data.list,PS_LIST_HEAD, false); + psMetadataItem *refItem = NULL; + while ((refItem = psListGetAndIncrement(iter))) { + if (refItem->type != PS_DATA_METADATA) { + psLogMsg ("psastro", PS_LOG_INFO, "GLINT.FILTER.TERM entry not metadata"); + return false; + } + char *refFilter = psMetadataLookupStr (&status, refItem->data.md, "FILTER"); + if (!status) { + // psLogMsg ("psastro", PS_LOG_INFO, "a PHOTCODE.DATA recipe folder is missing FILTER"); + continue; + } + if (strcmp(refFilter, filter)) continue; + + C_terms = psMetadataLookupPtr(&status, refItem->data.md, "C_TERMS"); + R_terms = psMetadataLookupPtr(&status, refItem->data.md, "R_TERMS"); + scale_factor = psMetadataLookupF32(&status, refItem->data.md, "FACTOR"); + double limit = psMetadataLookupF32(&status, refItem->data.md, "LIMIT"); + if (star->Mag > limit) { continue; } + } + // check + double glintPixelScale = psMetadataLookupF32(&status, glintItem->data.md, "GLINT.PIXEL.SCALE"); + // psMetadata *GLINT_PARAMETER_SET = psMetadataLookupPtr(&status, glintItem->data.md, "GLINT.FILTER.TERM"); + // psMetadata *GLINT_PARAMETERS = psMetadataLookupPtr(&status, GLINT_PARAMETER_SET, filter); + + + // double theta0 = atan2(star->FP->y,star->FP->x); + double theta0 = atan2(star->TP->y,star->TP->x); + double star_radius_deg = R_FPA * glintPixelScale; + + double C = C_terms->data.F32[0] + C_terms->data.F32[1] * star_radius_deg + C_terms->data.F32[2] * pow(star_radius_deg,2); + double R = R_terms->data.F32[0] + R_terms->data.F32[1] * star_radius_deg + R_terms->data.F32[2] * pow(star_radius_deg,2); + + double inner_edge_angle, outer_edge_angle; + if (star_radius_deg < 0.8831) { + outer_edge_angle = 0.0; + } + else if (star_radius_deg < 0.9368) { + outer_edge_angle = 46.2985 * sqrt(star_radius_deg - 0.8831); + } + else { + outer_edge_angle = -2.67 + 14.3 * star_radius_deg; + } + if (star_radius_deg < 0.964) { + inner_edge_angle = -243.866 * pow(star_radius_deg - 0.932,2) + 2.4636; + } + else { + inner_edge_angle = 66.2061 * sqrt(star_radius_deg - 0.9635); + } + + inner_edge_angle = inner_edge_angle * PS_RAD_DEG; + outer_edge_angle = outer_edge_angle * PS_RAD_DEG; + + // These define the ends of a single arc in arc-centric coordinates + double x1 = C - R * cos(inner_edge_angle); + double y1 = R * sin(inner_edge_angle); + double x2 = C - R * cos(outer_edge_angle); + double y2 = R * sin(outer_edge_angle); + + // Create the endpoints for the pair of arcs in device coordinates (in millimeters) + double theta1 = theta0 + M_PI / 2.0; + double x1s = (x1 * cos(theta1) - y1 * sin(theta1)) * scale_factor; + double y1s = (y1 * cos(theta1) + x1 * sin(theta1)) * scale_factor; + double x1e = (x2 * cos(theta1) - y2 * sin(theta1)) * scale_factor; + double y1e = (y2 * cos(theta1) + x2 * sin(theta1)) * scale_factor; + + double x2s = (x1 * cos(theta1) + y1 * sin(theta1)) * scale_factor; + double y2s = (-1.0 * y1 * cos(theta1) + x1 * sin(theta1)) * scale_factor; + double x2e = (x2 * cos(theta1) + y2 * sin(theta1)) * scale_factor; + double y2e = (-1.0 * y2 * cos(theta1) + x2 * sin(theta1)) * scale_factor; + + psVector *x_start = psVectorAlloc(4,PS_TYPE_F32); + psVector *y_start = psVectorAlloc(4,PS_TYPE_F32); + psVector *x_end = psVectorAlloc(4,PS_TYPE_F32); + psVector *y_end = psVectorAlloc(4,PS_TYPE_F32); + + // CZW: I know this looks like a typo, but trust me, it's not a typo. + psPlane *fp = psPlaneAlloc(); + psPlane *tp = psPlaneAlloc(); + + fprintf(stderr,"HSC GLINT: x1x2(%f %f) (%f %f) x1sx1e (%f %f %f %f) x2sx2e (%f %f %f %f)\n", + x1,y1,x2,y2, + x1s,y1s,x1e,y1e, + x2s,y2s,x2e,y2e); + + + tp->x = 13.5 * (y1s / 0.015 + 12.78); + tp->y = 13.5 * (x1s / -0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + x_start->data.F32[0] = fp->x; + y_start->data.F32[0] = fp->y; + x_start->data.F32[0] = tp->x; + y_start->data.F32[0] = tp->y; + + tp->x = 13.5 * (y1e / 0.015 + 12.78); + tp->y = 13.5 * (x1e / -0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + x_end->data.F32[0] = fp->x; + y_end->data.F32[0] = fp->y; + x_end->data.F32[0] = tp->x; + y_end->data.F32[0] = tp->y; + + x_start->data.F32[1] = x_end->data.F32[0]; + y_start->data.F32[1] = y_end->data.F32[0]; + x_end->data.F32[1] = x_start->data.F32[0]; + y_end->data.F32[1] = y_start->data.F32[0]; + + tp->x = 13.5 * (y2s / 0.015 + 12.78); + tp->y = 13.5 * (x2s / -0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + x_start->data.F32[2] = fp->x; + y_start->data.F32[2] = fp->y; + x_start->data.F32[2] = tp->x; + y_start->data.F32[2] = tp->y; + + tp->x = 13.5 * (y2e / 0.015 + 12.78); + tp->y = 13.5 * (x2e / -0.015 + 57.74); + psPlaneTransformApply(fp,fpa->fromTPA, tp); + x_end->data.F32[2] = fp->x; + y_end->data.F32[2] = fp->y; + x_end->data.F32[2] = tp->x; + y_end->data.F32[2] = tp->y; + + x_start->data.F32[3] = x_end->data.F32[2]; + y_start->data.F32[3] = y_end->data.F32[2]; + x_end->data.F32[3] = x_start->data.F32[2]; + y_end->data.F32[3] = y_start->data.F32[2]; + + psFree(fp); + psFree(tp); + + for (int glint_point = 0; glint_point < 4; glint_point++) { + for (int nChip = 0; nChip < fpa->chips->n; nChip++) { + pmChip *chip = fpa->chips->data[nChip]; + if (!chip) continue; + + if (!psastroFindChipInYrange (fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point])) { + continue; + } + if (!psastroFindChipInXrange (fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point])) { + continue; + } + + double xChip0, yChip0, xChip1, yChip1, chip_angle, glint_length; + psastroFPAtoChip (&xChip0, &yChip0, fpa, nChip, x_start->data.F32[glint_point], y_start->data.F32[glint_point]); + psastroFPAtoChip (&xChip1, &yChip1, fpa, nChip, x_end->data.F32[glint_point], y_end->data.F32[glint_point]); + + chip_angle = atan2(yChip1 - yChip0, xChip1 - xChip0); + glint_length = sqrt(pow(yChip1 - yChip0,2) + pow(xChip1 - xChip0,2)); + + // select the 0th readout of the 0th cell for this chip + if (!chip->cells) continue; + if (!chip->cells->n) continue; + pmCell *glintCell = chip->cells->data[0]; + if (!glintCell) continue; + if (!glintCell->readouts) continue; + if (!glintCell->readouts->n) continue; + pmReadout *glintReadout = glintCell->readouts->data[0]; + if (!glintReadout) continue; + + // save the glints on the readout->analysis metadata, creating if needed + psArray *glints = psMetadataLookupPtr (&status, glintReadout->analysis, "PSASTRO.GLINTS.HSC"); + if (glints == NULL) { + glints = psArrayAllocEmpty (100); + if (!psMetadataAdd (glintReadout->analysis, PS_LIST_TAIL, "PSASTRO.GLINTS.HSC", PS_DATA_ARRAY, "astrometry matches", glints)) { + psWarning("failure to add glints to readout"); + psFree (glints); + continue; + } + psFree (glints); + } + + fprintf (stderr, "glint %s : %d %f,%f to %f,%f (%f %f %f)\n", glintType, nChip, xChip0, yChip0, xChip1, yChip1, glint_length, glintWidth, chip_angle); + psTrace("psastro.masks",4,"HSC_GLINT: Star: %f %f Glint CR %f %f Parameters: %f %f %f Ends: %f %f -> %f %f Chip: %f %f -> %f %f @ %s %f\n", + star->FP->x,star->FP->y, + C,R, + inner_edge_angle,outer_edge_angle,theta0, + x_start->data.F32[glint_point],y_start->data.F32[glint_point],x_end->data.F32[glint_point],y_end->data.F32[glint_point], + xChip0,yChip0,xChip1,yChip1, + psMetadataLookupStr(&status,glintReadout->parent->parent->concepts,"CHIP.NAME"),chip_angle); + psVector *glint = psVectorAlloc(5,PS_TYPE_F32); + glint->data.F32[0] = xChip0; + glint->data.F32[1] = yChip0; + glint->data.F32[2] = glint_length * pixelScale / pixelScale; + glint->data.F32[3] = glintWidth; + glint->data.F32[4] = chip_angle; + + psArrayAdd (glints, 100, glint); + + psFree (glint); + } // End loop over chips + } // End loop over glint endpoints + + psFree(x_start); + psFree(x_end); + psFree(y_start); + psFree(y_end); + } // End HSC glint block + } } Index: branches/czw_branch/20170908/psastro/src/psastroMaskUpdates.c =================================================================== --- branches/czw_branch/20170908/psastro/src/psastroMaskUpdates.c (revision 40129) +++ branches/czw_branch/20170908/psastro/src/psastroMaskUpdates.c (revision 40179) @@ -129,5 +129,8 @@ bool REFSTAR_MASK_BLEED = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED"); - + bool REFSTAR_MASK_BLEED_ORIENTATION_X = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED_ORIENTATION_X"); + bool REFSTAR_MASK_BLEED_ORIENTATION_Y = psMetadataLookupBool (&status, recipe, "REFSTAR_MASK_BLEED_ORIENTATION_Y"); + + double REFSTAR_MASK_MAX_MAG = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_MAX_MAG"); double REFSTAR_MASK_SATSTAR_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_SATSTAR_MAG_SLOPE"); @@ -143,4 +146,9 @@ double REFSTAR_MASK_BLEED_MAG_SLOPE = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE"); + double REFSTAR_MASK_BLEED_MAG_MAX_X = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_MAX_X"); + double REFSTAR_MASK_BLEED_MAG_SLOPE_X = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE_X"); + double REFSTAR_MASK_BLEED_MAG_MAX_Y = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_MAX_Y"); + double REFSTAR_MASK_BLEED_MAG_SLOPE_Y = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_BLEED_MAG_SLOPE_Y"); + double REFSTAR_MASK_CROSSTALK_MAG_MAX = psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_MAX"); double REFSTAR_MASK_CROSSTALK_MAG_SLOPE= psMetadataLookupF32 (&status, recipe, "REFSTAR_MASK_CROSSTALK_MAG_SLOPE"); @@ -180,4 +188,6 @@ } + psTrace("psastro.masks",2,"Configuration loaded for masking."); + // recipe values are given in instrumental magnitudes // use the zero point and exposure time to convert to apparent mags: M_ap = M_inst + C_0 + 2.5*log(exptime) @@ -187,6 +197,12 @@ REFSTAR_MASK_SATSPIKE_MAG_MAX += MagOffset; REFSTAR_MASK_BLEED_MAG_MAX += MagOffset; + REFSTAR_MASK_BLEED_MAG_MAX_X += MagOffset; + REFSTAR_MASK_BLEED_MAG_MAX_Y += MagOffset; + GHOST_MAX_MAG += MagOffset; + psTrace("psastro.masks",2,"Magnitudes: max_mag: %f satstar: %f satspike: %f bleed: %f ghost %f\n", + REFSTAR_MASK_MAX_MAG,REFSTAR_MASK_SATSTAR_MAG_MAX,REFSTAR_MASK_SATSPIKE_MAG_MAX,REFSTAR_MASK_BLEED_MAG_MAX,GHOST_MAX_MAG); + // select the output mask image :: we mosaic to chip mosaic format pmFPAfile *outMask = psMetadataLookupPtr (NULL, config->files, "PSASTRO.OUTPUT.MASK"); @@ -287,4 +303,5 @@ if (!refstars) continue; + psTrace("psastro.masks",2,"In Readout loop."); // we need to generate the following masks regions: // 1) circle around the saturated stars (scaled by magnitude) @@ -301,5 +318,6 @@ } } - + psTrace("psastro.masks",4,"In refstar loop: %d/%ld %f %f\n", + i,refstars->n,ref->Mag,REFSTAR_MASK_MAX_MAG); if (ref->Mag > REFSTAR_MASK_MAX_MAG) continue; @@ -334,5 +352,7 @@ // The length should also be a function of the image background level - + psTrace("psastro.masks",4,"Masking: Radius: %f Theta: %f Length: %f Width: %f\n", + radius,Theta,spikeLength,spikeWidth); + psTrace("psastro.masks",4,"Also: %f %f %f %f\n",ref->chip->x,ref->chip->y,ROTANGLE,REFSTAR_MASK_SATSTAR_POS_ZERO); psastroMaskBox (readoutMask->mask, spikeMaskValue, ref->chip->x, ref->chip->y, spikeLength, spikeWidth, Theta); } @@ -349,15 +369,33 @@ float xCell = 0.0; float yCell = 0.0; - float xEnd = 0.0; - float yEnd = 0.0; + float xStart = 0.0; + float xEnd = 0.0; + float yStart = 0.0; + float yEnd = 0.0; + // float xEnd = 0.0; + // float yEnd = 0.0; + float width = 0.0; + float length = 0.0; // find coordinate of star on cell pmCellCoordsForChip (&xCell, &yCell, refCell, ref->chip->x, ref->chip->y); // find coordinate of end-point on chip - int ySize = psMetadataLookupS32(NULL, refCell->concepts, "CELL.YSIZE"); - pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell, ySize); - - float width = REFSTAR_MASK_BLEED_MAG_SLOPE*(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag); - psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) ref->chip->x-0.5*width, (int) ref->chip->y, (int) ref->chip->x+0.5*width+1, yEnd); + if (REFSTAR_MASK_BLEED_ORIENTATION_X) { + length = pow(10,REFSTAR_MASK_BLEED_MAG_SLOPE_X*(REFSTAR_MASK_BLEED_MAG_MAX_X - ref->Mag)); + width = REFSTAR_MASK_BLEED_MAG_SLOPE *(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag); + + pmChipCoordsForCell (&xStart, &yStart, refCell, xCell - length, yCell - 0.5 * width); + pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell + length, yCell + 0.5 * width); + psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) xStart, (int) yStart, (int) xEnd, (int) yEnd + 1); + } + if (REFSTAR_MASK_BLEED_ORIENTATION_Y) { + length = pow(10,REFSTAR_MASK_BLEED_MAG_SLOPE_Y*(REFSTAR_MASK_BLEED_MAG_MAX_Y - ref->Mag)); + width = REFSTAR_MASK_BLEED_MAG_SLOPE *(REFSTAR_MASK_BLEED_MAG_MAX - ref->Mag); + + pmChipCoordsForCell (&xStart, &yStart, refCell, xCell - 0.5 * width, yCell - length); + pmChipCoordsForCell (&xEnd, &yEnd, refCell, xCell + 0.5 * width, yCell + length); + psastroMaskRectangle (readoutMask->mask, spikeMaskValue, (int) xStart, (int) yStart, (int) xEnd + 1, (int) yEnd); + } + } } @@ -392,5 +430,21 @@ } } - + // Because it's very GPC1-specific, handle alternate glint types here + psArray *hsc_glints = psMetadataLookupPtr (&status, readout->analysis, "PSASTRO.GLINTS.HSC"); + if (hsc_glints) { + for (int i = 0; i < hsc_glints->n; i++) { + // Contains (start_x, start_y, length, width, theta) + psVector *glint = hsc_glints->data[i]; + + psTrace("psastro.masks", 4, "Masking glint %d on chip %s with params (%f,%f,%f,%f,%f)\n", + i, psMetadataLookupStr(&status,readout->parent->parent->concepts,"CHIP.NAME"), + glint->data.F32[0], glint->data.F32[1], + glint->data.F32[2], glint->data.F32[3], glint->data.F32[4]); + + psastroMaskBox (readoutMask->mask, glintMaskValue, glint->data.F32[0], glint->data.F32[1], + glint->data.F32[2], glint->data.F32[3], glint->data.F32[4]); + } + } + // Select the crosstalk artifact regions for this readout, and mask a circular region // corresponding to the source star