source: trunk/Mars/fact/processing/magicweather.C@ 20000

Last change on this file since 20000 was 18946, checked in by Daniela Dorner, 7 years ago
allow values +-45sec from the run
File size: 3.1 KB
Line 
1void magicweather(const char *fname, double beg=0, double end=100000)
2{
3 fits file(fname);
4
5 //file.PrintColumns();
6 //file.PrintKeys();
7
8 Double_t time;
9 Float_t temp, hum, pres, dew, wind, gust;
10 file.SetPtrAddress("Time", &time);
11 file.SetPtrAddress("T", &temp);
12 file.SetPtrAddress("H", &hum);
13 file.SetPtrAddress("P", &pres);
14 file.SetPtrAddress("T_dew", &dew);
15 file.SetPtrAddress("v", &wind);
16 file.SetPtrAddress("v_max", &gust);
17
18 UInt_t offset = file.GetUInt("MJDREF");
19 if (beg < 30000)
20 beg+=offset;
21 if (end < 30000)
22 end+=offset;
23
24 int cnt = 0;
25 double avg_t = 0;
26 double rms_t = 0;
27 double avg_h = 0;
28 double rms_h = 0;
29 double avg_p = 0;
30 double rms_p = 0;
31 double avg_d = 0;
32 double rms_d = 0;
33 double avg_w = 0;
34 double rms_w = 0;
35 double avg_g = 0;
36 double rms_g = 0;
37
38 double last_t = -1;
39 double last_h = -1;
40 double last_p = -1;
41 double last_d = -1;
42 double last_w = -1;
43 double last_g = -1;
44 double diff = -1;
45 double diff2 = -1;
46
47 while (file.GetNextRow())
48 {
49 time += offset;
50
51 if (time>end)
52 {
53 diff2 = time-end;
54 if (diff2<diff)
55 {
56 last_t = temp;
57 last_h = hum;
58 last_p = pres;
59 last_d = dew;
60 last_w = wind;
61 last_g = gust;
62 diff=diff2;
63 }
64 break;
65 }
66
67 if (time<beg)
68 {
69 last_t = temp;
70 last_h = hum;
71 last_p = pres;
72 last_d = dew;
73 last_w = wind;
74 last_g = gust;
75 diff = beg-time;
76 continue;
77 }
78
79 avg_t += temp;
80 rms_t += temp*temp;
81 avg_h += hum;
82 rms_h += hum*hum;
83 avg_p += pres;
84 rms_p += pres*pres;
85 avg_d += dew;
86 rms_d += dew*dew;
87 avg_w += wind;
88 rms_w += wind*wind;
89 avg_g += gust;
90 rms_g += gust*gust;
91 cnt ++;
92 }
93
94 if (cnt==0)
95 {
96 //only give output if values within 45 sec before/after are available
97 // weather reports come every ~30 sec
98 if (diff<0 || diff>45./24/3600)
99 return;
100
101 cout << "result " << last_t << " 0 " << last_h << " 0 "
102 << last_p << " 0 " << last_d << " 0 "
103 << last_w << " 0 " << last_g << " 0 " << endl;
104 return;
105 }
106
107 avg_t /= cnt;
108 rms_t /= cnt;
109 avg_h /= cnt;
110 rms_h /= cnt;
111 avg_p /= cnt;
112 rms_p /= cnt;
113 avg_d /= cnt;
114 rms_d /= cnt;
115 avg_w /= cnt;
116 rms_w /= cnt;
117 avg_g /= cnt;
118 rms_g /= cnt;
119
120 rms_t = sqrt(rms_t-avg_t*avg_t);
121 rms_h = sqrt(rms_h-avg_h*avg_h);
122 rms_p = sqrt(rms_p-avg_p*avg_p);
123 rms_d = sqrt(rms_d-avg_d*avg_d);
124 rms_w = sqrt(rms_w-avg_w*avg_w);
125 rms_g = sqrt(rms_g-avg_g*avg_g);
126
127 cout << "result " << avg_t << " " << rms_t << " " << avg_h << " " << rms_h << " "
128 << avg_p << " " << rms_p << " " << avg_d << " " << rms_d << " "
129 << avg_w << " " << rms_w << " " << avg_g << " " << rms_g << " " << endl;
130}
Note: See TracBrowser for help on using the repository browser.