1 #include <clocale>
2 #include <cstdio>
3 #include <cstdlib>
4 #include <fstream>
5 #include <sstream>
6 using namespace std;
7
8 #include "findFile.h"
9 #include "readOriginFile.h"
10 #include "xpUtil.h"
11
12 static double
interpolateCircular(double x0,double x1,const double xMax,const double interpFactor)13 interpolateCircular(double x0, double x1,
14 const double xMax, const double interpFactor)
15 {
16 const double larger = (x0 > x1 ? x0 : x1);
17 const double smaller = (x0 < x1 ? x0 : x1);
18
19 if ((larger - smaller) > (smaller - larger + xMax))
20 {
21 if (x0 < x1)
22 x0 += xMax;
23 else
24 x1 += xMax;
25 }
26
27 double returnVal = interpFactor * (x1 - x0) + x0;
28 if (returnVal > xMax)
29 returnVal -= xMax;
30
31 return(returnVal);
32 }
33
34 void
readOriginFile(string filename,vector<LBRPoint> & originVector)35 readOriginFile(string filename, vector<LBRPoint> &originVector)
36 {
37 if (!findFile(filename, "origin"))
38 {
39 ostringstream errStr;
40 errStr << "Can't open origin file " << filename << "\n";
41 xpExit(errStr.str(), __FILE__, __LINE__);
42 }
43
44 checkLocale(LC_NUMERIC, "C");
45
46 ifstream infile(filename.c_str());
47 char line[256];
48 while(infile.getline(line, 256))
49 {
50 if (line[0] == '#') continue;
51
52 long int yyyymmdd, hhmmss;
53 double r, lat, lon, localTime;
54
55 int numRead = sscanf(line, "%ld.%ld %lf %lf %lf %lf",
56 &yyyymmdd, &hhmmss,
57 &r, &lat, &lon, &localTime);
58 int yyyymm = yyyymmdd / 100;
59 int year = yyyymm/100;
60 int month = abs(yyyymm - year * 100);
61 int day = abs((int) yyyymmdd - yyyymm * 100);
62
63 int hhmm = hhmmss / 100;
64 int hour = hhmm / 100;
65 int min = hhmm - hour * 100;
66 int sec = hhmmss - hhmm * 100;
67
68 const double julian_day = toJulian(year, month, day, hour, min, sec);
69
70 // If only the date has been specified, set range to 0
71 if (numRead < 3) r = 0;
72
73 // If localTime isn't specified, set it to -1
74 if (numRead < 6) localTime = -1;
75
76 lat *= deg_to_rad;
77 lon *= deg_to_rad;
78
79 LBRPoint p = { julian_day, r, lat, lon, localTime };
80
81 originVector.push_back(p);
82 }
83 checkLocale(LC_NUMERIC, "");
84 infile.close();
85 }
86
87 void
readDynamicOrigin(string filename,LBRPoint & originPoint)88 readDynamicOrigin(string filename, LBRPoint &originPoint)
89 {
90 if (!findFile(filename, "origin"))
91 {
92 ostringstream errStr;
93 errStr << "Can't open origin file " << filename << "\n";
94 xpExit(errStr.str(), __FILE__, __LINE__);
95 }
96
97 checkLocale(LC_NUMERIC, "C");
98 ifstream infile(filename.c_str());
99 char line[256];
100 while(infile.getline(line, 256))
101 {
102 if (line[0] == '#') continue;
103
104 long int yyyymmdd, hhmmss; // not used
105 double r, lat, lon, localTime = -1;
106 sscanf(line, "%ld.%ld %lf %lf %lf %lf", &yyyymmdd, &hhmmss,
107 &r, &lat, &lon, &localTime);
108 lat *= deg_to_rad;
109 lon *= deg_to_rad;
110
111 originPoint.time = 0;
112 originPoint.radius = r;
113 originPoint.latitude = lat;
114 originPoint.longitude = lon;
115 originPoint.localTime = localTime;
116 }
117 checkLocale(LC_NUMERIC, "");
118 infile.close();
119 }
120
121 void
interpolateOriginFile(const double julianDay,const vector<LBRPoint> & originVector,double & rad,double & lat,double & lon,double & localTime)122 interpolateOriginFile(const double julianDay,
123 const vector<LBRPoint> &originVector,
124 double &rad,
125 double &lat, double &lon,
126 double &localTime)
127 {
128 const int numPoints = originVector.size();
129
130 if (julianDay < originVector[0].time)
131 {
132 rad = originVector[0].radius;
133 lat = originVector[0].latitude;
134 lon = originVector[0].longitude;
135 localTime = originVector[0].localTime;
136 }
137 else if (julianDay > originVector[numPoints-1].time)
138 {
139 rad = originVector[numPoints-1].radius;
140 lat = originVector[numPoints-1].latitude;
141 lon = originVector[numPoints-1].longitude;
142 localTime = originVector[numPoints-1].localTime;
143 }
144 else
145 {
146 double interpFactor = 0;
147 for (int i = 1; i < numPoints; i++)
148 {
149 if (julianDay < originVector[i].time)
150 {
151 const double jd0 = originVector[i-1].time;
152 const double jd1 = originVector[i].time;
153
154 interpFactor = (julianDay - jd0)/(jd1 - jd0);
155
156 rad = (interpFactor
157 * (originVector[i].radius - originVector[i-1].radius)
158 + originVector[i-1].radius);
159 lat = (interpFactor
160 * (originVector[i].latitude
161 - originVector[i-1].latitude)
162 + originVector[i-1].latitude);
163 lon = interpolateCircular(originVector[i-1].longitude,
164 originVector[i].longitude,
165 TWO_PI,
166 interpFactor);
167 localTime = interpolateCircular(originVector[i-1].localTime,
168 originVector[i].localTime,
169 24,
170 interpFactor);
171
172 break;
173 }
174 }
175 }
176 }
177