AstroLibrary

1
//------------------------------------------------------------------------------
2
//
3
// File:    Coco.cpp
4
//
5
// Coordinate transformations
6
//
7
//    The original Fortran/Pascal/C++ codes is described in
8
//
9
//  - Shampine, Gordon: "Computer solution of Ordinary Differential Equations",
10
//    Freeman and Comp., San Francisco (1975)
11
//  - Montenbruck O., Pfleger T.: Astronomie Mit Dem Personal Computer.
12
//    Springer-Verlag Berlin Heidelberg (1989).
13
//  - Montenbruck O., Pfleger T.: Astronomy on the Personal Computer.
14
//    Springer-Verlag Berlin Heidelberg (2000).
15
//------------------------------------------------------------------------------
16

17
#include <iomanip>
18
//#include <stdio.h>      // Header file for standard Input/Output
19

20
#include "APC_Const.h"
21
#include "APC_Math.h"
22
#include "APC_PrecNut.h"
23
#include "APC_Spheric.h"
24
#include "APC_Sun.h"
25
#include "APC_Time.h"
26
#include "APC_VecMat3D.h"
27

28
#ifdef __GNUC__   // GNU C++ Adaptations
29
#include <ctype.h>
30
#include "GNU_iomanip.h"
31
#endif
32

33
#pragma argsused
34

35
using namespace std;
36

37
//
38
// Types
39
//
40
enum enOrigin   { Heliocentric, Geocentric };  // Origin of coordinates
41
enum enRefSys   { Ecliptic, Equator };         // Reference system
42

43
//
44
// Definition of class "Position"
45
//
46
class Position
47
{
48
  public:
49
	void Input(); // Query user for parameters
50
	void SetOrigin(enOrigin Origin);
51
	void SetRefSys(enRefSys RefSys);
52
	void SetEquinox(double T_Equinox);
53
	void Print();
54
  private:
55
	Vec3D      m_R;          // Coordinate vector
56
	enOrigin   m_Origin;     // Origin of coordinate system
57
	enRefSys   m_RefSys;     // Reference system
58
	double     m_TEquinox;   // Equinox (cent. since J2000)
59
	double     m_MjdEpoch;   // Epoch (Modif. Julian Date)
60
};
61

62

63
//
64
// Data input
65
//
66
void Position::Input()
67
{
68
  //
69
  // Variables
70
  //
71
  char c;
72

73

74
  // Header
75
  cout << endl << "New input:" << endl << endl;
76

77

78
  // Query reference system
79
  while (true) {
80

81
    cout << "  Reference system (e=ecliptic,a=equator) ... ";
82
    cin >> c; cin.ignore(81,'\n'); c = tolower(c);
83

84
    if (c=='e') { m_RefSys=Ecliptic; break; };
85
    if (c=='a') { m_RefSys=Equator;  break; };
86
  }
87

88

89
  // Query coordinates
90
  while (true) {
91

92
    cout << "  Format (c=cartesian,p=polar)            ... ";
93
    cin >> c; cin.ignore(81,'\n'); c = tolower(c);
94

95
    if (c=='c') {
96
      double x,y,z;
97
      cout << "  Coordinates (x y z)                     ... ";
98
      cin >> x >> y >> z;
99
      m_R = Vec3D(x,y,z);
100
	  break;
101
    };
102

103
    if (c=='p') {
104

105
      if (m_RefSys==Ecliptic) {
106

107
        int     d,m;
108
        double  s, L,B,R;
109
        cout << "  Coordinates (L [o ' \"] B[o ' \"] R)     ... ";
110
        cin >> d >> m >> s; L=Rad*Ddd(d,m,s);
111
        cin >> d >> m >> s; B=Rad*Ddd(d,m,s);
112
        cin >> R; cin.ignore(81,'\n');
113
        m_R = Vec3D(Polar(L,B,R));
114
      }
115

116
      else {
117

118
        int     d,m;
119
		double  s, Ra,Dec,R;
120
        cout << "  Coordinates (Ra [h m s] Dec [o ' \"] R)  ... ";
121
        cin >> d >> m >> s; Ra=15.0*Rad*Ddd(d,m,s);
122
        cin >> d >> m >> s; Dec=Rad*Ddd(d,m,s);
123
        cin >> R; cin.ignore(81,'\n');
124
		m_R = Vec3D(Polar(Ra,Dec,R));
125
      };
126
	  break;
127
	};
128
  }
129

130

131
  // Query equinox
132
  double Year;
133

134
  cout << "  Equinox (yyyy.y)                        ... ";
135
  cin >> Year; cin.ignore(81,'\n');
136

137
  m_TEquinox = (Year-2000.0)/100.0;
138

139
  // Query origin
140
  while (true) {
141

142
    cout << "  Origin (h=heliocentric,g=geocentric)    ... ";
143
    cin >> c; cin.ignore(81,'\n'); c = tolower(c);
144

145
    if (c=='g') { m_Origin=Geocentric;   break; };
146
    if (c=='h') { m_Origin=Heliocentric; break; };
147
  }
148

149

150
  // Query epoch
151
  int    year,month,day;
152
  double Hour;
153

154
  cout << "  Epoch (yyyy mm dd hh.h)                 ... ";
155
  cin >> year >> month >> day >> Hour; cin.ignore(81,'\n');
156
  m_MjdEpoch = Mjd(year,month,day)+Hour/24.0;
157
  cout << endl; // Trailer
158
}
159

160

161
//
162
// Change of origin
163
//
164
void Position::SetOrigin(enOrigin Origin)
165
{
166
  //
167
  // Variables
168
  //
169
  double T_Epoch;
170
  Vec3D  R_Sun;
171

172
  if (Origin!=m_Origin) {
173
	// Geocentric coordinates of the Sun at epoch w.r.t.
174
    // given reference system and equinox
175
    T_Epoch = (m_MjdEpoch-MJD_J2000)/36525.0;
176
    if (m_RefSys==Ecliptic)
177
      R_Sun = PrecMatrix_Ecl(T_Epoch,m_TEquinox) * SunPos(T_Epoch);
178
    else
179
      R_Sun = Ecl2EquMatrix(m_TEquinox) *
180
              PrecMatrix_Ecl(T_Epoch,m_TEquinox) * SunPos(T_Epoch);
181
    // Change origin
182
	if (m_Origin==Heliocentric) {
183
      m_R += R_Sun; m_Origin = Geocentric;
184
    }
185
    else {
186
      m_R -= R_Sun; m_Origin = Heliocentric;
187
    };
188
  };
189
}
190

191
//
192
// Change of reference system
193
//
194
void Position::SetRefSys(enRefSys RefSys)
195
{
196
  if (RefSys!=m_RefSys) {
197
    if (m_RefSys==Equator) {
198
      m_R = Equ2EclMatrix(m_TEquinox) * m_R;  m_RefSys = Ecliptic;
199
	}
200
    else {
201
      m_R = Ecl2EquMatrix(m_TEquinox) * m_R;  m_RefSys = Equator;
202
    }
203
  };
204
}
205

206
//
207
// Change of equinox
208
//
209
void Position::SetEquinox(double T_Equinox)
210
{
211
  if (T_Equinox!=m_TEquinox) {
212
    if (m_RefSys==Equator)
213
      m_R = PrecMatrix_Equ(m_TEquinox,T_Equinox) * m_R;
214
    else
215
      m_R = PrecMatrix_Ecl(m_TEquinox,T_Equinox) * m_R;
216
    m_TEquinox = T_Equinox;
217
  };
218
}
219

220
//
221
// Output
222
//
223
void Position::Print()
224
{
225
  cout << endl;
226
  cout << (( m_Origin==Heliocentric )?  "Heliocentric " : "Geocentric ");
227
  cout << (( m_RefSys==Equator )?  "equatorial " : "ecliptic ");
228
  cout << "coordinates" << endl;
229

230
  cout << "(Equinox J" << fixed << setprecision(1)
231
       << 2000.0+m_TEquinox*100.0 << ",  ";
232

233
  cout << "Epoch " << DateTime(m_MjdEpoch,HHh) << ")" << endl;
234
  cout << endl;
235

236
  cout << "  (x,y,z) = " << fixed << setprecision(8) << setw(12) << m_R << endl;
237

238
  if ( m_RefSys==Equator ) {
239
    cout << "         h  m  s               o  '  \"" << endl;
240
    cout << "  Ra = " << setprecision(2) << setw(11)
241
         << Angle(Deg*m_R[phi]/15.0,DMMSSs);
242
    cout << "    Dec = " << setprecision(1) << showpos << setw(11)
243
         << Angle(Deg*m_R[theta],DMMSSs) << noshowpos;
244
  }
245
  else {
246
    cout << "         o  '  \"             o  '  \"" << endl;
247
    cout << "  L = " << setprecision(2) << setw(12)
248
         << Angle(Deg*m_R[phi],DMMSSs);
249
    cout << "    B = " << setprecision(1) << showpos << setw(11)
250
         << Angle(Deg*m_R[theta],DMMSSs) << noshowpos;
251
  }
252
  cout << "    R = " << setprecision(8) << setw(12)
253
       << m_R[r] << endl;
254
  cout << endl << endl;
255
}
256

257

258
//------------------------------------------------------------------------------
259
//
260
// Main program
261
//
262
//------------------------------------------------------------------------------
263
void main(int argc, _TCHAR* argv[]) {
264

265
  //
266
  // Variables
267
  //
268
  Position  Pos;
269
  char      c;
270
  bool      End = false;
271
  double    Year;
272

273
  // Header
274
  cout << endl
275
       << "                COCO: coordinate conversions       " << endl
276
       << "        (c) 1999 Oliver Montenbruck, Thomas Pfleger" << endl
277
       << endl;
278

279
  // Initialization
280
  Pos.Input();
281
  Pos.Print();
282

283

284
  // Command loop
285
  do {
286
	// Command input
287
	cout << "Enter command (?=Help) ... ";
288
	cin >> c; cin.ignore(81,'\n'); c=tolower(c);
289

290
	// Actions
291
	switch (c) {
292
	  case 'x':
293
		End = true;  break;
294
	  case 'a':
295
		Pos.SetRefSys(Equator); Pos.Print();  break;
296
	  case 'e':
297
		Pos.SetRefSys(Ecliptic); Pos.Print();  break;
298
	  case 'g':
299
		Pos.SetOrigin(Geocentric); Pos.Print();  break;
300
	  case 'h':
301
		Pos.SetOrigin(Heliocentric); Pos.Print();  break;
302
	  case 'n':
303
		Pos.Input(); Pos.Print();  break;
304
	  case 'p':
305
		cout << "New equinox (yyyy.y)   ... ";
306
		cin >> Year; cin.ignore(81,'\n');
307
		Pos.SetEquinox( (Year-2000.0)/100.0 );
308
		Pos.Print();
309
		break;
310
	  default:
311
		// Display help text
312
		cout << endl
313
			 << "Available commands" << endl
314
			 << "  a=equatorial, e=ecliptic, p=precession,     " << endl
315
			 << "  g=geocentric, h=heliocentric, n=new input,  " << endl
316
			 << "  x=exit                                      " << endl
317
			 << endl;
318
		break;
319
	}
320
  }
321
  while (!End);
322

323
  cout << endl;
324
}
325