~ubuntu-branches/ubuntu/oneiric/gkrellmoon/oneiric : revision 1

1

/*

2

* Calculates Moon's position according to Brown's Lunar Theory...

3

*

4

* (C) Mike Henderson <mghenderson@lanl.gov>.

5

*

6

* I've converted to glib types, removed unused variables, and piped the whole

7

* thing through indent.

8

*

9

* josh buhl <jbuhl@users.sourceforge.net>

10

*/

11

12

#include <math.h>

13

#include "Moon.h"

14

15

#define DegPerRad 57.29577951308232087680

16

#define RadPerDeg 0.01745329251994329576

17

18

gdouble TwoPi = 6.283185308;

19

gdouble ARC = 206264.81;

20

gdouble DLAM, DLAMS;

21

gdouble DS;

22

gdouble GAM1C;

23

gdouble SINPI;

24

gdouble N;

25

gdouble CO[14][5], SI[14][5];

26

27

28

static gdouble sine(gdouble phi)

29

{

30

return (sin(TwoPi * frac(phi)));

31

32

}

33

34

35

gdouble frac(gdouble x)

36

{

37

38

x -= (gint) x;

39

return ((x < 0) ? x + 1.0 : x);

40

41

}

42

43

44

static void

45

addthe(gdouble C1, gdouble S1, gdouble C2, gdouble S2, gdouble * C,

46

gdouble * S)

47

{

48

49

*C = C1 * C2 - S1 * S2;

50

*S = S1 * C2 + C1 * S2;

51

52

53

}

54

55

56

static void term(gint P, gint Q, gint R, gint S, gdouble * X, gdouble * Y)

57

{

58

59

gdouble XX, YY;

60

gint k, I[5];

61

62

I[1] = P, I[2] = Q, I[3] = R, I[4] = S, XX = 1.0, YY = 0.0;

63

for (k = 1; k <= 4; ++k) {

64

if (I[k] != 0.0) {

65

addthe(XX, YY, CO[6 + I[k]][k], SI[6 + I[k]][k],

66

&XX, &YY);

67

}

68

}

69

*X = XX;

70

*Y = YY;

71

72

}

73

74

75

static void

76

addsol(gdouble COEFFL, gdouble COEFFS, gdouble COEFFG, gdouble COEFFP,

77

gint P, gint Q, gint R, gint S)

78

{

79

80

gdouble X, Y;

81

82

term(P, Q, R, S, &X, &Y);

83

DLAM += COEFFL * Y;

84

DS += COEFFS * Y;

85

GAM1C += COEFFG * X;

86

SINPI += COEFFP * X;

87

88

89

}

90

91

92

93

static void addn(gdouble COEFFN, gint P, gint Q, gint R, gint S)

94

{

95

96

gdouble X, Y;

97

98

term(P, Q, R, S, &X, &Y);

99

N += COEFFN * Y;

100

101

102

}

103

104

105

gdouble

106

Moon(gdouble T, gdouble * LAMBDA, gdouble * BETA, gdouble * R,

107

gdouble * AGE)

108

{

109

110

gdouble T2;

111

gdouble S1, S2, S3, S4, S5, S6, S7;

112

gdouble DL0, DL, DD, DGAM, DLS, DF;

113

gdouble L0, L, LS, F, D;

114

gdouble ARG, FAC;

115

gint MAX, i, j;

116

gdouble S;

117

118

ARG = 0.0, FAC = 0.0, MAX = 0;

119

120

T2 = T * T;

121

DLAM = 0.0, DS = 0.0, GAM1C = 0.0;

122

SINPI = 3422.7000;

123

124

/*

125

* Long Periodic variations

126

*/

127

S1 = sine(0.19833 + 0.05611 * T);

128

S2 = sine(0.27869 + 0.04508 * T);

129

S3 = sine(0.16827 - 0.36903 * T);

130

S4 = sine(0.34734 - 5.37261 * T);

131

S5 = sine(0.10498 - 5.37899 * T);

132

S6 = sine(0.42681 - 0.41855 * T);

133

S7 = sine(0.14943 - 5.37511 * T);

134

DL0 =

135

0.84 * S1 + 0.31 * S2 + 14.27 * S3 + 7.26 * S4 + 0.28 * S5 +

136

0.24 * S6;

137

DL =

138

2.94 * S1 + 0.31 * S2 + 14.27 * S3 + 9.34 * S4 + 1.12 * S5 +

139

0.83 * S6;

140

DLS = -6.40 * S1 - 1.89 * S6;

141

DF =

142

0.21 * S1 + 0.31 * S2 + 14.27 * S3 - 88.70 * S4 - 15.30 * S5 +

143

0.24 * S6 - 1.86 * S7;

144

DD = DL0 - DLS;

145

DGAM = -3332e-9 * sine(0.59734 - 5.37261 * T)

146

- 539e-9 * sine(0.35498 - 5.37899 * T)

147

- 64e-9 * sine(0.39943 - 5.37511 * T);

148

149

150

151

L0 =

152

TwoPi * frac(0.60643382 + 1336.85522467 * T -

153

0.00000313 * T2) + DL0 / ARC;

154

L =

155

TwoPi * frac(0.37489701 + 1325.55240982 * T +

156

0.00002565 * T2) + DL / ARC;

157

LS =

158

TwoPi * frac(0.99312619 + 99.99735956 * T - 0.00000044 * T2) +

159

DLS / ARC;

160

F =

161

TwoPi * frac(0.25909118 + 1342.22782980 * T -

162

0.00000892 * T2) + DF / ARC;

163

D =

164

TwoPi * frac(0.82736186 + 1236.85308708 * T -

165

0.00000397 * T2) + DD / ARC;

166

167

for (i = 1; i <= 4; ++i) {

168

switch (i) {

169

case 1:

170

ARG = L, MAX = 4, FAC = 1.000002208;

171

break;

172

case 2:

173

ARG = LS, MAX = 3, FAC =

174

0.997504612 - 0.002495388 * T;

175

break;

176

case 3:

177

ARG = F, MAX = 4, FAC =

178

1.000002708 + 139.978 * DGAM;

179

break;

180

case 4:

181

ARG = D, MAX = 6, FAC = 1.0;

182

break;

183

}

184

185

CO[6 + 0][i] = 1.0, CO[6 + 1][i] = cos(ARG) * FAC;

186

SI[6 + 0][i] = 0.0, SI[6 + 1][i] = sin(ARG) * FAC;

187

for (j = 2; j <= MAX; ++j)

188

addthe(CO[6 + j - 1][i], SI[6 + j - 1][i],

189

CO[6 + 1][i], SI[6 + 1][i], &CO[6 + j][i],

190

&SI[6 + j][i]);

191

for (j = 1; j <= MAX; ++j) {

192

CO[6 - j][i] = CO[6 + j][i];

193

SI[6 - j][i] = -SI[6 + j][i];

194

}

195

196

197

}

198

199

200

201

/*

202

* Solar1

203

*/

204

addsol(13.902, 14.06, -0.001, 0.2607, 0, 0, 0, 4);

205

addsol(0.403, -4.01, +0.394, 0.0023, 0, 0, 0, 3);

206

addsol(2369.912, 2373.36, +0.601, 28.2333, 0, 0, 0, 2);

207

addsol(-125.154, -112.79, -0.725, -0.9781, 0, 0, 0, 1);

208

addsol(1.979, 6.98, -0.445, 0.0433, 1, 0, 0, 4);

209

addsol(191.953, 192.72, +0.029, 3.0861, 1, 0, 0, 2);

210

addsol(-8.466, -13.51, +0.455, -0.1093, 1, 0, 0, 1);

211

addsol(22639.500, 22609.07, +0.079, 186.5398, 1, 0, 0, 0);

212

addsol(18.609, 3.59, -0.094, 0.0118, 1, 0, 0, -1);

213

addsol(-4586.465, -4578.13, -0.077, 34.3117, 1, 0, 0, -2);

214

addsol(+3.215, 5.44, +0.192, -0.0386, 1, 0, 0, -3);

215

addsol(-38.428, -38.64, +0.001, 0.6008, 1, 0, 0, -4);

216

addsol(-0.393, -1.43, -0.092, 0.0086, 1, 0, 0, -6);

217

addsol(-0.289, -1.59, +0.123, -0.0053, 0, 1, 0, 4);

218

addsol(-24.420, -25.10, +0.040, -0.3000, 0, 1, 0, 2);

219

addsol(18.023, 17.93, +0.007, 0.1494, 0, 1, 0, 1);

220

addsol(-668.146, -126.98, -1.302, -0.3997, 0, 1, 0, 0);

221

addsol(0.560, 0.32, -0.001, -0.0037, 0, 1, 0, -1);

222

addsol(-165.145, -165.06, +0.054, 1.9178, 0, 1, 0, -2);

223

addsol(-1.877, -6.46, -0.416, 0.0339, 0, 1, 0, -4);

224

addsol(0.213, 1.02, -0.074, 0.0054, 2, 0, 0, 4);

225

addsol(14.387, 14.78, -0.017, 0.2833, 2, 0, 0, 2);

226

addsol(-0.586, -1.20, +0.054, -0.0100, 2, 0, 0, 1);

227

addsol(769.016, 767.96, +0.107, 10.1657, 2, 0, 0, 0);

228

addsol(+1.750, 2.01, -0.018, 0.0155, 2, 0, 0, -1);

229

addsol(-211.656, -152.53, +5.679, -0.3039, 2, 0, 0, -2);

230

addsol(+1.225, 0.91, -0.030, -0.0088, 2, 0, 0, -3);

231

addsol(-30.773, -34.07, -0.308, 0.3722, 2, 0, 0, -4);

232

addsol(-0.570, -1.40, -0.074, 0.0109, 2, 0, 0, -6);

233

addsol(-2.921, -11.75, +0.787, -0.0484, 1, 1, 0, 2);

234

addsol(+1.267, 1.52, -0.022, 0.0164, 1, 1, 0, 1);

235

addsol(-109.673, -115.18, +0.461, -0.9490, 1, 1, 0, 0);

236

addsol(-205.962, -182.36, +2.056, +1.4437, 1, 1, 0, -2);

237

addsol(0.233, 0.36, 0.012, -0.0025, 1, 1, 0, -3);

238

addsol(-4.391, -9.66, -0.471, 0.0673, 1, 1, 0, -4);

239

240

241

/*

242

* Solar2

243

*/

244

addsol(0.283, 1.53, -0.111, +0.0060, 1, -1, 0, +4);

245

addsol(14.577, 31.70, -1.540, +0.2302, 1, -1, 0, 2);

246

addsol(147.687, 138.76, +0.679, +1.1528, 1, -1, 0, 0);

247

addsol(-1.089, 0.55, +0.021, 0.0, 1, -1, 0, -1);

248

addsol(28.475, 23.59, -0.443, -0.2257, 1, -1, 0, -2);

249

addsol(-0.276, -0.38, -0.006, -0.0036, 1, -1, 0, -3);

250

addsol(0.636, 2.27, +0.146, -0.0102, 1, -1, 0, -4);

251

addsol(-0.189, -1.68, +0.131, -0.0028, 0, 2, 0, 2);

252

addsol(-7.486, -0.66, -0.037, -0.0086, 0, 2, 0, 0);

253

addsol(-8.096, -16.35, -0.740, 0.0918, 0, 2, 0, -2);

254

addsol(-5.741, -0.04, 0.0, -0.0009, 0, 0, 2, 2);

255

addsol(0.255, 0.0, 0.0, 0.0, 0, 0, 2, 1);

256

addsol(-411.608, -0.20, 0.0, -0.0124, 0, 0, 2, 0);

257

addsol(0.584, 0.84, 0.0, +0.0071, 0, 0, 2, -1);

258

addsol(-55.173, -52.14, 0.0, -0.1052, 0, 0, 2, -2);

259

addsol(0.254, 0.25, 0.0, -0.0017, 0, 0, 2, -3);

260

addsol(+0.025, -1.67, 0.0, +0.0031, 0, 0, 2, -4);

261

addsol(1.060, 2.96, -0.166, 0.0243, 3, 0, 0, +2);

262

addsol(36.124, 50.64, -1.300, 0.6215, 3, 0, 0, 0);

263

addsol(-13.193, -16.40, +0.258, -0.1187, 3, 0, 0, -2);

264

addsol(-1.187, -0.74, +0.042, 0.0074, 3, 0, 0, -4);

265

addsol(-0.293, -0.31, -0.002, 0.0046, 3, 0, 0, -6);

266

addsol(-0.290, -1.45, +0.116, -0.0051, 2, 1, 0, 2);

267

addsol(-7.649, -10.56, +0.259, -0.1038, 2, 1, 0, 0);

268

addsol(-8.627, -7.59, +0.078, -0.0192, 2, 1, 0, -2);

269

addsol(-2.740, -2.54, +0.022, 0.0324, 2, 1, 0, -4);

270

addsol(1.181, 3.32, -0.212, 0.0213, 2, -1, 0, +2);

271

addsol(9.703, 11.67, -0.151, 0.1268, 2, -1, 0, 0);

272

addsol(-0.352, -0.37, +0.001, -0.0028, 2, -1, 0, -1);

273

addsol(-2.494, -1.17, -0.003, -0.0017, 2, -1, 0, -2);

274

addsol(0.360, 0.20, -0.012, -0.0043, 2, -1, 0, -4);

275

addsol(-1.167, -1.25, +0.008, -0.0106, 1, 2, 0, 0);

276

addsol(-7.412, -6.12, +0.117, 0.0484, 1, 2, 0, -2);

277

addsol(-0.311, -0.65, -0.032, 0.0044, 1, 2, 0, -4);

278

addsol(+0.757, 1.82, -0.105, 0.0112, 1, -2, 0, 2);

279

addsol(+2.580, 2.32, +0.027, 0.0196, 1, -2, 0, 0);

280

addsol(+2.533, 2.40, -0.014, -0.0212, 1, -2, 0, -2);

281

addsol(-0.344, -0.57, -0.025, +0.0036, 0, 3, 0, -2);

282

addsol(-0.992, -0.02, 0.0, 0.0, 1, 0, 2, 2);

283

addsol(-45.099, -0.02, 0.0, -0.0010, 1, 0, 2, 0);

284

addsol(-0.179, -9.52, 0.0, -0.0833, 1, 0, 2, -2);

285

addsol(-0.301, -0.33, 0.0, 0.0014, 1, 0, 2, -4);

286

addsol(-6.382, -3.37, 0.0, -0.0481, 1, 0, -2, 2);

287

addsol(39.528, 85.13, 0.0, -0.7136, 1, 0, -2, 0);

288

addsol(9.366, 0.71, 0.0, -0.0112, 1, 0, -2, -2);

289

addsol(0.202, 0.02, 0.0, 0.0, 1, 0, -2, -4);

290

291

/*

292

* Solar3

293

*/

294

addsol(0.415, 0.10, 0.0, 0.0013, 0, 1, 2, 0);

295

addsol(-2.152, -2.26, 0.0, -0.0066, 0, 1, 2, -2);

296

addsol(-1.440, -1.30, 0.0, +0.0014, 0, 1, -2, 2);

297

addsol(0.384, -0.04, 0.0, 0.0, 0, 1, -2, -2);

298

addsol(+1.938, +3.60, -0.145, +0.0401, 4, 0, 0, 0);

299

addsol(-0.952, -1.58, +0.052, -0.0130, 4, 0, 0, -2);

300

addsol(-0.551, -0.94, +0.032, -0.0097, 3, 1, 0, 0);

301

addsol(-0.482, -0.57, +0.005, -0.0045, 3, 1, 0, -2);

302

addsol(0.681, 0.96, -0.026, 0.0115, 3, -1, 0, 0);

303

addsol(-0.297, -0.27, 0.002, -0.0009, 2, 2, 0, -2);

304

addsol(0.254, +0.21, -0.003, 0.0, 2, -2, 0, -2);

305

addsol(-0.250, -0.22, 0.004, 0.0014, 1, 3, 0, -2);

306

addsol(-3.996, 0.0, 0.0, +0.0004, 2, 0, 2, 0);

307

addsol(0.557, -0.75, 0.0, -0.0090, 2, 0, 2, -2);

308

addsol(-0.459, -0.38, 0.0, -0.0053, 2, 0, -2, 2);

309

addsol(-1.298, 0.74, 0.0, +0.0004, 2, 0, -2, 0);

310

addsol(0.538, 1.14, 0.0, -0.0141, 2, 0, -2, -2);

311

addsol(0.263, 0.02, 0.0, 0.0, 1, 1, 2, 0);

312

addsol(0.426, +0.07, 0.0, -0.0006, 1, 1, -2, -2);

313

addsol(-0.304, +0.03, 0.0, +0.0003, 1, -1, 2, 0);

314

addsol(-0.372, -0.19, 0.0, -0.0027, 1, -1, -2, 2);

315

addsol(+0.418, 0.0, 0.0, 0.0, 0, 0, 4, 0);

316

addsol(-0.330, -0.04, 0.0, 0.0, 3, 0, 2, 0);

317

318

N = 0.0;

319

addn(-526.069, 0, 0, 1, -2);

320

addn(-3.352, 0, 0, 1, -4);

321

addn(+44.297, +1, 0, 1, -2);

322

addn(-6.000, +1, 0, 1, -4);

323

addn(+20.599, -1, 0, 1, 0);

324

addn(-30.598, -1, 0, 1, -2);

325

addn(-24.649, -2, 0, 1, 0);

326

addn(-2.000, -2, 0, 1, -2);

327

addn(-22.571, 0, +1, 1, -2);

328

addn(+10.985, 0, -1, 1, -2);

329

330

DLAM +=

331

+0.82 * sine(0.7736 - 62.5512 * T) + 0.31 * sine(0.0466 -

332

125.1025 * T)

333

+ 0.35 * sine(0.5785 - 25.1042 * T) + 0.66 * sine(0.4591 +

334

1335.8075 *

335

T) +

336

0.64 * sine(0.3130 - 91.5680 * T) + 1.14 * sine(0.1480 +

337

1331.2898 * T)

338

+ 0.21 * sine(0.5918 + 1056.5859 * T) + 0.44 * sine(0.5784 +

339

1322.8595 *

340

T) +

341

0.24 * sine(0.2275 - 5.7374 * T) + 0.28 * sine(0.2965 +

342

2.6929 * T)

343

+ 0.33 * sine(0.3132 + 6.3368 * T);

344

345

346

*LAMBDA = 360.0 * frac((L0 + DLAM / ARC) / TwoPi);

347

348

S = F + DS / ARC;

349

FAC = 1.000002708 + 139.978 * DGAM;

350

*BETA =

351

(FAC * (18518.511 + 1.189 + GAM1C) * sin(S) -

352

6.24 * sin(3 * S) + N) / 3600.0;

353

354

SINPI *= 0.999953253;

355

*R = ARC / SINPI;

356

357

358

DLAMS = 6893.0 * sin(LS) + 72.0 * sin(2.0 * LS);

359

360

*AGE = 29.530589 * frac((D + (DLAM - DLAMS) / ARC) / TwoPi);

361

/*

362

printf("Diff = %f\n", 360.0*frac((D+(DLAM-DLAMS)/ARC)/TwoPi));

363

*/

364

365

/*

366

* Return the phase.

367

*/

368

/*

369

return( 0.5*(1.0 - cos(D+(DLAM-DLAMS)/ARC)) );

370

*/

371

return (*AGE / 29.530589);

372

373

374

375

}

376

377

#define R 0.61803399

378

#define C 0.38196601

379

380

gdouble NewMoon(gdouble ax, gdouble bx, gdouble cx)

381

{

382

383

gdouble f1, f2, x0, x1, x2, x3, Moon();

384

gdouble L, B, Rad, AGE, tol = 1e-7;

385

386

x0 = ax;

387

x3 = cx;

388

if (fabs(cx - bx) > fabs(bx - ax)) {

389

x1 = bx;

390

x2 = bx + C * (cx - bx);

391

} else {

392

x2 = bx;

393

x1 = bx - C * (bx - ax);

394

}

395

f1 = Moon(x1, &L, &B, &Rad, &AGE);

396

f2 = Moon(x2, &L, &B, &Rad, &AGE);

397

while (fabs(x3 - x0) > tol * (fabs(x1) + fabs(x2))) {

398

if (f2 < f1) {

399

x0 = x1;

400

x1 = x2;

401

x2 = R * x1 + C * x3;

402

f1 = f2;

403

f2 = Moon(x2, &L, &B, &Rad, &AGE);

404

} else {

405

x3 = x2;

406

x2 = x1;

407

x1 = R * x2 + C * x0;

408

f2 = f1;

409

f1 = Moon(x1, &L, &B, &Rad, &AGE);

410

}

411

}

412

if (f1 < f2) {

413

return (x1);

414

} else {

415

return (x2);

416

}

417

}

418

419

420

421

422

423

/*

424

* MINI_MOON: low precision lunar coordinates (approx. 5'/1')

425

* T : time in Julian centuries since J2000

426

* ( T=(JD-2451545)/36525 )

427

* RA : right ascension (in h; equinox of date)

428

* DEC: declination (in deg; equinox of date)

429

*

430

*/

431

gint MiniMoon(gdouble T, gdouble * RA, gdouble * DEC)

432

{

433

434

gdouble L0, L, LS, F, D, H, S, N, DL, CB, L_MOON, B_MOON, V, W, X,

435

Y, Z, RHO;

436

gdouble frac(), cosEPS, sinEPS, P2, ARC;

437

438

439

cosEPS = 0.91748;

440

sinEPS = 0.39778;

441

P2 = 6.283185307;

442

ARC = 206264.8062;

443

444

445

/*

446

* mean elements of lunar orbit

447

*/

448

L0 = frac(0.606433 + 1336.855225 * T); /* mean longitude Moon (in rev) */

449

L = P2 * frac(0.374897 + 1325.552410 * T); /* mean anomaly of the Moon */

450

LS = P2 * frac(0.993133 + 99.997361 * T); /* mean anomaly of the Sun */

451

D = P2 * frac(0.827361 + 1236.853086 * T); /* diff. longitude Moon-Sun */

452

F = P2 * frac(0.259086 + 1342.227825 * T); /* mean argument of latitude */

453

DL =

454

+22640.0 * sin(L) - 4586.0 * sin(L - 2.0 * D) +

455

2370.0 * sin(2.0 * D) + 769.0 * sin(2.0 * L)

456

- 668.0 * sin(LS) - 412.0 * sin(2.0 * F) -

457

212.0 * sin(2.0 * L - 2.0 * D) - 206.0 * sin(L + LS - 2.0 * D)

458

+ 192.0 * sin(L + 2.0 * D) - 165.0 * sin(LS - 2.0 * D) -

459

125.0 * sin(D) - 110.0 * sin(L + LS)

460

+ 148.0 * sin(L - LS) - 55.0 * sin(2.0 * F - 2.0 * D);

461

S = F + (DL + 412.0 * sin(2.0 * F) + 541.0 * sin(LS)) / ARC;

462

H = F - 2.0 * D;

463

N =

464

-526.0 * sin(H) + 44.0 * sin(L + H) - 31.0 * sin(-L + H) -

465

23.0 * sin(LS + H)

466

+ 11.0 * sin(-LS + H) - 25.0 * sin(-2.0 * L + F) +

467

21.0 * sin(-L + F);

468

L_MOON = P2 * frac(L0 + DL / 1296e3); /* in rad */

469

B_MOON = (18520.0 * sin(S) + N) / ARC; /* in rad */

470

471

/* equatorial coordinates */

472

CB = cos(B_MOON);

473

X = CB * cos(L_MOON);

474

V = CB * sin(L_MOON);

475

W = sin(B_MOON);

476

Y = cosEPS * V - sinEPS * W;

477

Z = sinEPS * V + cosEPS * W;

478

RHO = sqrt(1.0 - Z * Z);

479

*DEC = (360.0 / P2) * atan2(Z, RHO);

480

*RA = (48.0 / P2) * atan2(Y, X + RHO);

481

if (*RA < 0.0)

482

*RA += 24.0;

483

484

return (0);

485

486

487

}