You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
 
 
 
 
 

333 lines
6.7 KiB

  1. /*
  2. * Atpdec
  3. * Copyright (c) 2003 by Thierry Leconte (F4DWV)
  4. *
  5. * $Id$
  6. *
  7. * This library is free software; you can redistribute it and/or modify
  8. * it under the terms of the GNU Library General Public License as
  9. * published by the Free Software Foundation; either version 2 of
  10. * the License, or (at your option) any later version.
  11. *
  12. * This program is distributed in the hope that it will be useful,
  13. * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15. * GNU Library General Public License for more details.
  16. *
  17. * You should have received a copy of the GNU Library General Public
  18. * License along with this library; if not, write to the Free Software
  19. * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20. *
  21. */
  22. #include <stdio.h>
  23. #include <string.h>
  24. #include <sndfile.h>
  25. #include <math.h>
  26. #define REGORDER 3
  27. typedef struct
  28. {
  29. double cf[REGORDER + 1];
  30. } rgparam;
  31. static void
  32. rgcomp (double x[16], rgparam * rgpr)
  33. {
  34. /*{ 0.106,0.215,0.324,0.433,0.542,0.652,0.78,0.87 ,0.0 }; */
  35. const double y[9] =
  36. { 31.07, 63.02, 94.96, 126.9, 158.86, 191.1, 228.62, 255.0, 0.0 };
  37. extern void polyreg (const int m, const int n, const double x[],
  38. const double y[], double c[]);
  39. polyreg (REGORDER, 9, x, y, rgpr->cf);
  40. }
  41. static double
  42. rgcal (float x, rgparam * rgpr)
  43. {
  44. double y, p;
  45. int i;
  46. for (i = 0, y = 0.0, p = 1.0; i < REGORDER + 1; i++) {
  47. y += rgpr->cf[i] * p;
  48. p = p * x;
  49. }
  50. return (y);
  51. }
  52. static double tele[16];
  53. static double Cs;
  54. static nbtele;
  55. int
  56. Calibrate (float **prow, int nrow, int offset)
  57. {
  58. double teleline[3000];
  59. double wedge[16];
  60. rgparam regr[30];
  61. int n;
  62. int k;
  63. int mtelestart, telestart;
  64. int channel = -1;
  65. float max;
  66. printf ("Calibration ");
  67. fflush (stdout);
  68. /* build telemetry values lines */
  69. for (n = 0; n < nrow; n++) {
  70. int i;
  71. teleline[n] = 0.0;
  72. for (i = 3; i < 43; i++) {
  73. teleline[n] += prow[n][i + offset + 956];
  74. }
  75. teleline[n] /= 40.0;
  76. }
  77. if (nrow < 192) {
  78. fprintf (stderr, " impossible, not enought row\n");
  79. return (0);
  80. }
  81. /* find telemetry start in the 2nd third */
  82. max = 0.0;
  83. mtelestart = 0;
  84. for (n = nrow / 3 - 64; n < 2 * nrow / 3 - 64; n++) {
  85. float df;
  86. df =
  87. (teleline[n - 4] + teleline[n - 3] + teleline[n - 2] +
  88. teleline[n - 1]) / (teleline[n] + teleline[n + 1] + teleline[n + 2] +
  89. teleline[n + 3]);
  90. if (df > max) {
  91. mtelestart = n;
  92. max = df;
  93. }
  94. }
  95. mtelestart -= 64;
  96. telestart = mtelestart % 128;
  97. if (mtelestart < 0 || nrow < telestart + 128) {
  98. fprintf (stderr, " impossible, not enought row\n");
  99. return (0);
  100. }
  101. /* compute wedges and regression */
  102. for (n = telestart, k = 0; n < nrow - 128; n += 128, k++) {
  103. int j;
  104. for (j = 0; j < 16; j++) {
  105. int i;
  106. wedge[j] = 0.0;
  107. for (i = 1; i < 7; i++)
  108. wedge[j] += teleline[n + j * 8 + i];
  109. wedge[j] /= 6;
  110. }
  111. rgcomp (wedge, &(regr[k]));
  112. if (k == nrow / 256) {
  113. int i, l;
  114. /* telemetry calibration */
  115. for (j = 0; j < 16; j++) {
  116. tele[j] = rgcal (wedge[j], &(regr[k]));
  117. }
  118. /* channel ID */
  119. for (j = 0, max = 10000.0, channel = -1; j < 6; j++) {
  120. float df;
  121. df = wedge[15] - wedge[j];
  122. df = df * df;
  123. if (df < max) {
  124. channel = j;
  125. max = df;
  126. }
  127. }
  128. /* Cs computation */
  129. for (Cs = 0.0, i = 0, j = n; j < n + 128; j++) {
  130. double csline;
  131. for (csline = 0.0, l = 3; l < 43; l++)
  132. csline += prow[n][l + offset];
  133. csline /= 40.0;
  134. if (csline > 50.0) {
  135. Cs += csline;
  136. i++;
  137. }
  138. }
  139. Cs /= i;
  140. Cs = rgcal (Cs, &(regr[k]));
  141. }
  142. }
  143. nbtele = k;
  144. /* calibrate */
  145. for (n = 0; n < nrow; n++) {
  146. float *pixelv;
  147. int i, c;
  148. pixelv = prow[n];
  149. for (i = 0; i < 1001; i++) {
  150. float pv;
  151. int k, kof;
  152. pv = pixelv[i + offset];
  153. k = (n - telestart) / 128;
  154. if (k >= nbtele)
  155. k = nbtele - 1;
  156. kof = (n - telestart) % 128;
  157. if (kof < 64) {
  158. if (k < 1) {
  159. pv = rgcal (pv, &(regr[k]));
  160. }
  161. else {
  162. pv = rgcal (pv, &(regr[k])) * (64 + kof) / 128.0 +
  163. rgcal (pv, &(regr[k - 1])) * (64 - kof) / 128.0;
  164. }
  165. }
  166. else {
  167. if ((k + 1) >= nbtele) {
  168. pv = rgcal (pv, &(regr[k]));
  169. }
  170. else {
  171. pv = rgcal (pv, &(regr[k])) * (192 - kof) / 128.0 +
  172. rgcal (pv, &(regr[k + 1])) * (kof - 64) / 128.0;
  173. }
  174. }
  175. if (pv > 255.0)
  176. pv = 255.0;
  177. if (pv < 0.0)
  178. pv = 0.0;
  179. pixelv[i + offset] = pv;
  180. }
  181. }
  182. printf ("Done\n");
  183. return (channel);
  184. }
  185. /* ------------------------------temperature calibration -----------------------*/
  186. extern int satnum;
  187. const struct
  188. {
  189. float d[4][3];
  190. struct
  191. {
  192. float vc, A, B;
  193. } rad[3];
  194. struct
  195. {
  196. float Ns;
  197. float b[3];
  198. } cor[3];
  199. } satcal[2] =
  200. #include "satcal.h"
  201. typedef struct
  202. {
  203. double Nbb;
  204. double Cs;
  205. double Cb;
  206. int ch;
  207. } tempparam;
  208. /* temperature compensation for IR channel */
  209. static void
  210. tempcomp (double t[16], int ch, tempparam * tpr)
  211. {
  212. double Tbb, T[4];
  213. double C;
  214. double r;
  215. int n;
  216. tpr->ch = ch - 3;
  217. /* compute equivalent T black body */
  218. for (n = 0; n < 4; n++) {
  219. float d0, d1, d2;
  220. C = t[9 + n] * 4.0;
  221. d0 = satcal[satnum].d[n][0];
  222. d1 = satcal[satnum].d[n][1];
  223. d2 = satcal[satnum].d[n][2];
  224. T[n] = d0;
  225. T[n] += d1 * C;
  226. C = C * C;
  227. T[n] += d2 * C;
  228. }
  229. Tbb = (T[0] + T[1] + T[2] + T[3]) / 4.0;
  230. Tbb = satcal[satnum].rad[tpr->ch].A + satcal[satnum].rad[tpr->ch].B * Tbb;
  231. /* compute radiance Black body */
  232. C = satcal[satnum].rad[tpr->ch].vc;
  233. tpr->Nbb = c1 * C * C * C / (exp (c2 * C / Tbb) - 1.0);
  234. /* store Count Blackbody and space */
  235. tpr->Cs = Cs * 4.0;
  236. tpr->Cb = t[14] * 4.0;
  237. }
  238. static double
  239. tempcal (float Ce, tempparam * rgpr)
  240. {
  241. double Nl, Nc, Ns, Ne;
  242. double T, vc;
  243. Ns = satcal[satnum].cor[rgpr->ch].Ns;
  244. Nl = Ns + (rgpr->Nbb - Ns) * (rgpr->Cs - Ce * 4.0) / (rgpr->Cs - rgpr->Cb);
  245. Nc = satcal[satnum].cor[rgpr->ch].b[0] +
  246. satcal[satnum].cor[rgpr->ch].b[1] * Nl +
  247. satcal[satnum].cor[rgpr->ch].b[2] * Nl * Nl;
  248. Ne = Nl + Nc;
  249. vc = satcal[satnum].rad[rgpr->ch].vc;
  250. T = c2 * vc / log (c1 * vc * vc * vc / Ne + 1.0);
  251. T = (T - satcal[satnum].rad[rgpr->ch].A) / satcal[satnum].rad[rgpr->ch].B;
  252. /* rescale to range 0-255 for +40-60 °C */
  253. T = (-T + 273.15 + 40) / 100.0 * 255.0;
  254. return (T);
  255. }
  256. int
  257. Temperature (float **prow, int nrow, int channel, int offset)
  258. {
  259. tempparam temp;
  260. int n;
  261. printf ("Temperature ");
  262. tempcomp (tele, channel, &temp);
  263. for (n = 0; n < nrow; n++) {
  264. float *pixelv;
  265. int i, c;
  266. pixelv = prow[n];
  267. for (i = 0; i < 1001; i++) {
  268. float pv;
  269. pv = tempcal (pixelv[i + offset], &temp);
  270. if (pv > 255.0)
  271. pv = 255.0;
  272. if (pv < 0.0)
  273. pv = 0.0;
  274. pixelv[i + offset] = pv;
  275. }
  276. }
  277. printf ("Done\n");
  278. }