1 /********************************************************************
2 sweephgen4.c
3
4 Create ephemeris file type 4 ep4_
5 a fast precoomputed ephemeris used in some Astrodienst applications.
6
7 options: -fYYY (start) file number, required option
8 -nNN number of files to be created, default 1
9 -v verbose: print differences (default: no)
10 -t test by reading
11
12
13 File format:
14 1000 blocks of xxx bytes
15 File names: ep4_243, ep4_244
16 corresponding to the absolute julian day number
17
18 *********************************************************************/
19 /* Copyright (C) 1997 - 2021 Astrodienst AG, Switzerland. All rights reserved.
20
21 License conditions
22 ------------------
23
24 This file is part of Swiss Ephemeris.
25
26 Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
27 or distributor accepts any responsibility for the consequences of using it,
28 or for whether it serves any particular purpose or works at all, unless he
29 or she says so in writing.
30
31 Swiss Ephemeris is made available by its authors under a dual licensing
32 system. The software developer, who uses any part of Swiss Ephemeris
33 in his or her software, must choose between one of the two license models,
34 which are
35 a) GNU Affero General Public License (AGPL)
36 b) Swiss Ephemeris Professional License
37
38 The choice must be made before the software developer distributes software
39 containing parts of Swiss Ephemeris to others, and before any public
40 service using the developed software is activated.
41
42 If the developer choses the AGPL software license, he or she must fulfill
43 the conditions of that license, which includes the obligation to place his
44 or her whole software project under the AGPL or a compatible license.
45 See https://www.gnu.org/licenses/agpl-3.0.html
46
47 If the developer choses the Swiss Ephemeris Professional license,
48 he must follow the instructions as found in http://www.astro.com/swisseph/
49 and purchase the Swiss Ephemeris Professional Edition from Astrodienst
50 and sign the corresponding license contract.
51
52 The License grants you the right to use, copy, modify and redistribute
53 Swiss Ephemeris, but only under certain conditions described in the License.
54 Among other things, the License requires that the copyright notices and
55 this notice be preserved on all copies.
56
57 Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
58
59 The authors of Swiss Ephemeris have no control or influence over any of
60 the derived works, i.e. over software or services created by other
61 programmers which use Swiss Ephemeris functions.
62
63 The names of the authors or of the copyright holder (Astrodienst) must not
64 be used for promoting any software, product or service which uses or contains
65 the Swiss Ephemeris. This copyright notice is the ONLY place where the
66 names of the authors can legally appear, except in cases where they have
67 given special permission in writing.
68
69 The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
70 for promoting such software, products or services.
71 */
72
73
74 # include "swephexp.h"
75 # include "sweephe4.h"
76
77 # define EPHR_NPL (PLACALC_CHIRON + 1)
78
79 char *arg0;
80 int32 max_dd[EP_CALC_N]; /* remember maximum of second dfifferences */
81 double max_err[EP_CALC_N]; /* remember maximum error */
82 AS_BOOL verbose = FALSE;
83 char errtext[AS_MAXCH];
84
split(w,m,min,sec)85 int split(w, m, min, sec)
86 int32 w; /* position in seconds/m */
87 int m; /* factor for seconds */
88 short *min, /* storage for degrees and minutes */
89 *sec; /* storage for seconds * m */
90 {
91 if (w >= 0) {
92 *sec = w % (60 * m);
93 *min = w / (60 * m);
94 } else {
95 *sec = -(-w % (60 * m));
96 *min = -(-w / (60 * m));
97 }
98 return OK;
99 }
100
101
102 /*************************************************************
103 Pack positions of 10 days and write to file
104 The longitude is packed with second differences in such a way,
105 that the accumulating rounding erros do not exceed half of
106 the last stored digit, i.e. 0.05" moon, 0.005" other planets
107 **************************************************************/
eph4_pack(int32 jd,double (* l)[NDB],double ecliptic[],double nutation[])108 int eph4_pack (int32 jd, double (*l)[NDB], double ecliptic[],
109 double nutation[])
110 {
111 int i, p,ps;
112 int32 d1, d2, dd, d_ret, w0, w_ret;
113 double err;
114 struct ep4 e;
115 e.j_10000 = jd / 10000.0;
116 e.j_rest = jd - 10000.0 * e.j_10000;
117 w0 = swe_d2l( ecliptic[0] * DEG);
118 split( w0, 100, &e.ecl0m, &e.ecl0s );
119 for (i = 1; i < NDB; i++)
120 e.ecld1[i-1] = swe_d2l(ecliptic[i] * DEG - w0);
121 for (i = 0; i < NDB; i++)
122 e.nuts[i] = swe_d2l( nutation[i] * DEG ); /* int32 casted into short */
123 for (p = PLACALC_SUN; p <= PLACALC_CHIRON ; p++) {
124 ps = p;
125 w0 = swe_d2l( l[ps][0] * DEG);
126 d1 = swe_d2l( l[ps][1] * DEG - w0);
127 if (d1 >= DEG180)
128 d1 -= DEG360;
129 else if (d1 <= -DEG180)
130 d1 += DEG360;
131 split(w0, 100, &e.elo[p].p0m, &e.elo[p].p0s);
132 split(d1, 100, &e.elo[p].pd1m, &e.elo[p].pd1s);
133 d_ret = d1; /* recalculated diff */
134 w_ret = w0 + d_ret; /* recalculated position */
135 for (i = 2; i < NDB; i++) {
136 d2 = swe_d2l( l[ps][i] * DEG - w_ret);
137 if (d2 >= DEG180)
138 d2 -= DEG360;
139 else if (d2 <= -DEG180)
140 d2 += DEG360;
141 dd = d2 - d_ret; /* second difference */
142 if (p == PLACALC_MOON || p == PLACALC_MERCURY)
143 dd = swe_d2l(dd / 10.0); /* moon only 0.1" */
144 if (verbose && abs(dd) > abs(max_dd[ps]))
145 max_dd[ps] = dd;
146 e.elo[p].pd2[i-2] = dd;
147 if (p == PLACALC_MOON || p == PLACALC_MERCURY)
148 d_ret += e.elo[p].pd2[i-2] * 10L;
149 else
150 d_ret += e.elo[p].pd2[i-2];
151 w_ret += d_ret;
152 if (verbose) {
153 err = swe_difdeg2n(w_ret/360000.0, l[ps][i]); /* error */
154 if (fabs(err) > fabs(max_err[ps]))
155 max_err[ps] = err;
156 }
157 }
158 } /* for p */
159 #ifdef INTEL_BYTE_ORDER
160 shortreorder((UCHAR *) &e, sizeof(struct ep4));
161 #endif
162 fwrite (&e, sizeof(struct ep4), 1, ephfp);
163 return (OK);
164 }
165
166
167 /*************************************/
degstr(t)168 char *degstr (t)
169 double t;
170 {
171 static char a[20]; /* must survive call */
172 double min, sec;
173 int ideg, imin;
174 char sign = ' ';
175 if ( t < 0) sign = '-';
176 t = fabs (t);
177 ideg = (int) floor (t);
178 min = ( t - ideg ) * 60.0;
179 imin = (int) floor(min);
180 sec = ( min - imin ) * 60.0;
181 sprintf (a, "%c%3d %2d'%5.2f\"", sign, ideg, imin, sec);
182 return (a);
183 } /* degstr */
184
185 /********************************************************/
eph_test()186 int eph_test()
187 {
188 char cal;
189 int p, jday, jmon, jyear;
190 double al, jd;
191 centisec *cp;
192 while (TRUE) {
193 printf ("date ?");
194 if (scanf ("%d%d%d", &jday,&jmon,&jyear) < 1) exit(1);
195 if (jyear < 1600)
196 cal = 'j';
197 else
198 cal = 'g';
199 swe_date_conversion (jyear, jmon, jday, 0, cal, &jd);
200 if ((cp = ephread(jd, 0,0, errtext)) == NULL) {
201 fprintf (stderr,"%s: %s", arg0, errtext);
202 exit (1);
203 }
204 printf ("ephgen test d=%12.1f dmy %d.%d.%d", jd, jday, jmon, jyear);
205 if (cal == 'g')
206 printf (" greg");
207 else
208 printf (" julian");
209 printf ("\n\tecliptic %s ", degstr(cp[EP_ECL_INDEX]*CS2DEG));
210 printf ("nutation %s\n", degstr(cp[EP_NUT_INDEX] * CS2DEG));
211 for (p = 0; p <= PLACALC_CHIRON; p++) {
212 al = cp[p] * CS2DEG;
213 printf ("%2d%18s\n", p, degstr(al));
214 }
215 }
216 } /* end ephtest */
217
main(int argc,char ** argv)218 int main(int argc, char **argv)
219 {
220 int day, i, n, p;
221 char serr[AS_MAXCH];
222 double l[EPHR_NPL][NDB], ecliptic[NDB], nutation[NDB];
223 double jd0, jd;
224 double x[6];
225 int32 jlong;
226 int file;
227 int nfiles = 1;
228 int fnr = -10000;
229 int32 iflagret;
230 arg0 = argv[0];
231 for (i = 1; i < argc; i++) {
232 if (strncmp(argv[i], "-f", 2) == 0) {
233 fnr = atoi (argv[i] + 2);
234 if (fnr < -20 || fnr > 300) {
235 printf("file number out of range -20 ... 300");
236 exit (1);
237 }
238 }
239 if (strncmp(argv[i], "-n", 2) == 0) {
240 nfiles = atoi (argv[i] + 2);
241 }
242 if (strncmp(argv[i], "-v", 2) == 0) {
243 verbose = TRUE;
244 }
245 if (strncmp(argv[i], "-t", 2) == 0) {
246 eph_test();
247 exit(0);
248 }
249 }
250 if (fnr == -10000) {
251 fprintf(stderr,"missing file number -fNNN\n");
252 exit(1);
253 }
254 for (file = fnr; file < fnr + nfiles; file++) {
255 if (file > fnr) printf ("\n");
256 printf ("file = %d\n", file);
257 jd0 = EP4_NDAYS * file + 0.5;
258 jlong = floor(jd0);
259 if (eph4_posit (jlong, TRUE, errtext) != OK) {
260 fprintf (stderr,"%s: %s", arg0, errtext);
261 exit(1);
262 }
263 for (n = 0; n < EP4_NDAYS; n += NDB, jd0 += NDB) {
264 if (n % 500 == 0) {
265 if ( n > 0 && verbose) {
266 printf ("\ndd");
267 for (p = 0; p < 11; p++) {
268 printf("%6d ",max_dd[p]);
269 max_dd[p] = 0;
270 }
271 printf("\ner");
272 for (p = 0; p < 11; p++) {
273 printf("%6.3f ",max_err[p] * 3600);
274 max_err[p] = 0;
275 }
276 }
277 printf ("\n%d ", n);
278 } else {
279 printf (".");
280 }
281 fflush( stdout );
282 for (day = 0; day < NDB; day++) { /* compute positions for 10 days */
283 jd = jd0 + day;
284 for (p = PLACALC_SUN; p <= EP_CALC_N; p++) {
285 if ((iflagret = swe_calc(jd, ephe_plac2swe(p), 0, x, serr)) == ERR) {
286 swe_close();
287 printf("error in swe_calc() %s\n", serr);
288 exit (1);
289 }
290 l[p][day] = x[0];
291 }
292 if ((iflagret = swe_calc(jd, SE_ECL_NUT, 0, x, serr)) == ERR) {
293 swe_close();
294 printf("error in swe_calc() %s\n", serr);
295 exit (1);
296 }
297 ecliptic[day] = x[0];
298 nutation[day] = x[2];
299 }
300 jlong = floor(jd0);
301 eph4_pack (jlong, l, ecliptic, nutation);
302 }
303 putchar('\n');
304 fclose (ephfp);
305 ephfp = NULL;
306 } /* for file */
307 swe_close();
308 return(0);
309 } /* end main */
310