1 //

2 // Created by alex on 11/25/23.

3 //

4 #include "sun2000.h"

5 #include "nav.h"

6 double dpsi, eps;

7 int nutime = -9999;

8 void sun2000(int iyr, int iday, double sec, float sun[3], float* rs) {

9  const static int imon = 1;

10  const double static xk = 0.0056932; // !Constant of aberration

11  // Compute floating point days since Jan 1.5, 2000

12  // Note that the Julian day starts at noon on the specified date

13  double t = jd(iyr, imon, iday) - 2451545.0e0 + (sec - 43200.e0) / 86400.e0;

14  double xls, gs, xlm, omega;

15  // extern double dpsi, eps;

16  // extern int nutime;

17  // Compute solar ephemeris parameters

18  ephparms(t, &xls, &gs, &xlm, &omega);

19  // Check if need to compute nutation corrections for this day

20  int nt = (int)t;

21  if (nt != nutime) {

22  nutime = nt;

23  nutate(t, xls, gs, xlm, omega, &dpsi, &eps);

24  }

25 

26  // c Compute planet mean anomalies

27  // c Venus Mean Anomaly

28  double g2 = 50.40828e0 + 1.60213022e0 * t;

29  g2 = dmod(g2, 360.e0);

30 

31  // c Mars Mean Anomaly

32  double g4 = 19.38816e0 + 0.52402078e0 * t;

33  g4 = dmod(g4, 360.e0);

34 

35  // c Jupiter Mean Anomaly

36  double g5 = 20.35116e0 + 0.08309121e0 * t;

37  g5 = dmod(g5, 360.e0);

38 

39  // Compute solar distance(AU)

40  *rs = 1.00014e0 - 0.01671e0 * cos(gs / radeg) - 0.00014e0 * cos(2.0e0 * gs / radeg);

41 

42  // c Compute Geometric Solar Longitude

43  double dls = (6893.0e0 - 4.6543463e-4 * t) * sin(gs / radeg) + 72.0e0 * sin(2.0e0 * gs / radeg) -

44  7.0e0 * cos((gs - g5) / radeg) + 6.0e0 * sin((xlm - xls) / radeg) +

45  5.0e0 * sin((4.0e0 * gs - 8.0e0 * g4 + 3.0e0 * g5) / radeg) -

46  5.0e0 * cos((2.0e0 * gs - 2.0e0 * g2) / radeg) - 4.0e0 * sin((gs - g2) / radeg) +

47  4.0e0 * cos((4.0e0 * gs - 8.0e0 * g4 + 3.0e0 * g5) / radeg) +

48  3.0e0 * sin((2.0e0 * gs - 2.0e0 * g2) / radeg) - 3.0e0 * sin(g5 / radeg) -

49  3.0e0 * sin((2.0e0 * gs - 2.0e0 * g5) / radeg);

50  double xlsg = xls + dls / 3600.e0;

51  // c Compute Apparent Solar Longitude; includes corrections for nutation

52  // in longitude and velocity aberration

53  double xlsa = xlsg + dpsi - xk / (*rs);

54  // Compute unit Sun vector

55  sun[0] = cos(xlsa / radeg);

56  sun[1] = sin(xlsa / radeg) * cos(eps / radeg);

57  sun[2] = sin(xlsa / radeg) * sin(eps / radeg);

58 }