source: branches/removing_cpp11_features/fact/processing/currents.C@ 18407

Last change on this file since 18407 was 17100, checked in by Daniela Dorner, 12 years ago
added 3 new variables to be filled to the DB
File size: 3.5 KB
Line 
1
2
3void currents(const char *fname, double beg=0, double end=100000)
4{
5 fits file(fname);
6
7 //file.PrintColumns();
8 //file.PrintKeys();
9
10 Double_t time;
11 Double_t jd;
12
13 Float_t med, dev;
14 file.SetPtrAddress("Time", &time);
15 file.SetPtrAddress("I_med", &med);
16 file.SetPtrAddress("I_dev", &dev);
17
18 UInt_t offset = file.GetUInt("MJDREF");
19 if (beg < 30000)
20 beg+=offset;
21 if (end < 30000)
22 end+=offset;
23
24 double med_beg = 0;
25 double med_end = 0;
26 double med_avg = 0;
27 double med_rms = 0;
28 double dev_avg = 0;
29 double dev_rms = 0;
30 int cnt_beg = 0;
31 int cnt_end = 0;
32
33 double med_last = -1;
34 double dev_last = -1;
35 double diff = -1;
36
37 double calc_dev_avg = 0;
38 double calc_dev_rel = 0;
39
40 TGraph g;
41
42 while (file.GetNextRow())
43 {
44 time += offset;
45 jd = time + 2400000.5;
46
47 if (time>end)
48 break;
49
50 if (time<beg)
51 {
52 med_last = med;
53 dev_last = dev;
54 diff = beg-time;
55 continue;
56 }
57
58 med_avg += med;
59 med_rms += med*med;
60 dev_avg += dev;
61 dev_rms += dev*dev;
62
63 if (time<beg+15./24/3600)
64 {
65 med_beg += med;
66 cnt_beg++;
67 }
68 if (time>end-15./24/3600)
69 {
70 med_end += med;
71 cnt_end++;
72 }
73
74 //calculate the moon
75 Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);
76 // get local position of moon
77 Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);
78
79 //calculate the light condition
80 double disk = Nova::GetLunarDisk(jd);
81
82 // Current prediction
83 //double cang = sin(angle);
84 const double calt = sin(hrzm.alt*TMath::DegToRad());
85
86 // semi-major axis
87 const double lc = /*cang==0 ? 0 :*/ calt*pow(disk, 2.2)*pow(Nova::GetLunarEarthDist(jd)/384400, -2);///sqrt(cang);
88 //const double lc = cang==0 ? -1 : calt*pow(disk, 2.2)*pos(dist, -2);
89 const double Ical = lc>0 ? 4+103*lc : 4;
90
91 g.SetPoint(g.GetN(), time, med);
92
93 calc_dev_avg += med - Ical;
94 calc_dev_rel += med/Ical;
95
96 }
97
98 if (g.GetN()==0)
99 {
100 if (diff<5./24/3600)
101 return;
102
103 // return the last report before the run
104 // (if not older than 5 minutes)
105 cout << "result " << med_last << " 0 " << dev_last << " 0 0 0 0" << endl;
106 return;
107 }
108
109 Double_t a0, a1;
110 Int_t ifail;
111 g.LeastSquareLinearFit(0, a0, a1, ifail);
112
113 double fluct_abs = 0;
114 double fluct_rel = 0;
115 if (ifail==0)
116 {
117 for (int i=0; i<g.GetN(); i++)
118 g.GetY()[i] -= a0+a1*g.GetX()[i];
119 fluct_abs = g.GetRMS(2);
120
121 for (int i=0; i<g.GetN(); i++)
122 g.GetY()[i] /= a0+a1*g.GetX()[i];
123 fluct_rel = g.GetRMS(2);
124 }
125
126 med_beg = cnt_beg>0 ? med_beg/cnt_beg : -1;
127 med_end = cnt_end>0 ? med_end/cnt_end : -1;
128
129 med_avg /= g.GetN();
130 med_rms /= g.GetN();
131 med_rms = sqrt(med_rms-med_avg*med_avg);
132
133 dev_avg /= g.GetN();
134 dev_rms /= g.GetN();
135 dev_rms = sqrt(dev_rms-dev_avg*dev_avg);
136 calc_dev_avg /= g.GetN();
137
138 calc_dev_rel -= g.GetN();
139 calc_dev_rel /= g.GetN();
140
141
142 cout << "result " << med_avg << " " << med_rms << " " << dev_avg << " " << dev_rms << " " << med_beg << " " << med_end << " " << calc_dev_avg << " " << calc_dev_rel << " " << fluct_abs << " " << fluct_rel << endl;
143}
Note: See TracBrowser for help on using the repository browser.