source: Sophya/trunk/SophyaExt/XephemAstroLib/bdl.c@ 3444

Last change on this file since 3444 was 3111, checked in by cmv, 19 years ago

mise en conformite xephem 3.7.2 cmv 22/11/2006
File size: 5.3 KB
Line 
1/* crack natural satellite files from BDL */
2
3#include <stdio.h>
4#include <stdlib.h>
5#include <string.h>
6#include <math.h>
7#include <unistd.h>
8
9#include "astro.h"
10#include "bdl.h"
11
12typedef enum {I, F, NL} ScanType;
13#define SCANFLD(f,w,vp) if(readField(fp,f,w,(void *)vp,ynot)<0) return (-1)
14
15static int readField (FILE *fp, ScanType f, int w, void *ptr, char ynot[]);
16static int readRec (FILE *fp, double *t0, double cmx[], double cfx[],
17 double cmy[], double cfy[], double cmz[], double cfz[], char ynot[]);
18
19/* given a sequencial text file in BDL natural satellite ephemeris format and a
20 * JD, find the x/y/z positions of each satellite. store in the given arrays,
21 * assumed to have one entry per moon. values are planetocentric, +x east, +y
22 * north, +z away from earth, all in au. corrected for light time.
23 * return the number of satellites or -1 and reason in ymot[].
24 * files obtained from ftp://ftp.bdl.fr/pub/misc/satxyz.
25 */
26int
27read_bdl (FILE *fp, double jd, double *xp, double *yp, double *zp, char ynot[])
28{
29 int npla;
30 int nsat;
31 int idn[8];
32 double freq[8];
33 double delt[8];
34 double djj;
35 double cmx[6], cfx[4], cmy[6], cfy[4], cmz[6], cfz[4];
36 int ienrf;
37 int jan;
38 int reclen;
39 long os0;
40 double t0;
41 int i;
42
43 /* read header line */
44 SCANFLD (I, 2, &npla);
45 SCANFLD (I, 2, &nsat);
46 for (i = 0; i < nsat; i++)
47 SCANFLD (I, 5, &idn[i]);
48 for (i = 0; i < nsat; i++)
49 SCANFLD (F, 8, &freq[i]);
50 for (i = 0; i < nsat; i++)
51 SCANFLD (F, 5, &delt[i]);
52 SCANFLD (I, 5, &ienrf);
53 SCANFLD (F, 15, &djj);
54 SCANFLD (I, 5, &jan);
55 SCANFLD (NL, 0, NULL);
56
57 /* record position of first record */
58 os0 = ftell (fp);
59
60 /* read first record to get length */
61 reclen = readRec (fp, &t0, cmx, cfx, cmy, cfy, cmz, cfz, ynot);
62 if (reclen < 0)
63 return (-1);
64
65 /* compute location of each satellite */
66 for (i = 0; i < nsat; i++) {
67 int id = (int)floor((jd-djj)/delt[i]) + idn[i] - 2;
68 long os = os0 + id*reclen;
69 double t1, anu, tau, tau2, at;
70 double tbx, tby, tbz;
71
72 if (fseek (fp, os, SEEK_SET) < 0) {
73 sprintf (ynot, "Seek error to %ld for rec %d", os, id);
74 return (-1);
75 }
76
77 if (readRec (fp, &t0, cmx, cfx, cmy, cfy, cmz, cfz, ynot) < 0)
78 return (-1);
79
80 t1 = floor(t0) + 0.5;
81 anu = freq[i];
82 tau = jd - t1;
83 tau2 = tau * tau;
84 at = tau*anu;
85
86 tbx = cmx[0]+cmx[1]*tau+cmx[2]*sin(at+cfx[0])
87 +cmx[3]*tau*sin(at+cfx[1])
88 +cmx[4]*tau2*sin(at+cfx[2])
89 +cmx[5]*sin(2*at+cfx[3]);
90 tby = cmy[0]+cmy[1]*tau+cmy[2]*sin(at+cfy[0])
91 +cmy[3]*tau*sin(at+cfy[1])
92 +cmy[4]*tau2*sin(at+cfy[2])
93 +cmy[5]*sin(2*at+cfy[3]);
94 tbz = cmz[0]+cmz[1]*tau+cmz[2]*sin(at+cfz[0])
95 +cmz[3]*tau*sin(at+cfz[1])
96 +cmz[4]*tau2*sin(at+cfz[2])
97 +cmz[5]*sin(2*at+cfz[3]);
98
99 xp[i] = tbx*1000./149597870.;
100 yp[i] = tby*1000./149597870.;
101 zp[i] = tbz*1000./149597870.;
102 }
103
104 return (nsat);
105}
106
107/* read one field.
108 * return 0 if ok else -1
109 * N.B. this is enforce width, without skipping leading blanks.
110 */
111static int
112readField (FILE *fp, ScanType f, int width, void *ptr, char ynot[])
113{
114 char buf[128];
115 char *bp;
116
117 if (buf && width != (int)fread (buf, 1, width, fp)) {
118 if (ferror(fp)) strcpy (ynot, "IO error");
119 else if (feof(fp)) strcpy (ynot, "unexpected EOF");
120 else strcpy (ynot, "short file");
121 return (-1);
122 }
123
124 buf[width] = '\0';
125 switch (f) {
126 case I:
127 *(int *)ptr = atoi (buf);
128 break;
129 case F:
130 bp = strchr (buf, 'D');
131 if (bp)
132 *bp = 'e';
133 *(double *)ptr = atof (buf);
134 break;
135 case NL:
136 fgets (buf, sizeof(buf), fp);
137 break;
138 default:
139 sprintf (ynot, "Bug! format = %d", f);
140 return (-1);
141 }
142 return (0);
143}
144
145/* read one satellite record.
146 * return number of chars read else -1.
147 */
148static int
149readRec (FILE *fp, double *t0, double cmx[], double cfx[], double cmy[],
150double cfy[], double cmz[], double cfz[], char ynot[])
151{
152
153 long pos0, pos1;
154 int isat, idx;
155 int ldat1, ldat2;
156 int i;
157
158 pos0 = ftell (fp);
159
160 SCANFLD (I, 1, &isat);
161 SCANFLD (I, 5, &idx);
162 SCANFLD (I, 8, &ldat1);
163 SCANFLD (I, 8, &ldat2);
164 SCANFLD (F, 9, t0);
165 for (i = 0; i < 6; i++)
166 SCANFLD (F, 17, &cmx[i]);
167 for (i = 0; i < 4; i++)
168 SCANFLD (F, 17, &cfx[i]);
169 for (i = 0; i < 6; i++)
170 SCANFLD (F, 17, &cmy[i]);
171 for (i = 0; i < 4; i++)
172 SCANFLD (F, 17, &cfy[i]);
173 for (i = 0; i < 6; i++)
174 SCANFLD (F, 17, &cmz[i]);
175 for (i = 0; i < 4; i++)
176 SCANFLD (F, 17, &cfz[i]);
177 SCANFLD (NL, 0, NULL);
178
179 pos1 = ftell (fp);
180
181 return (pos1 - pos0);
182}
183
184
185#ifdef TEST_IT
186/* stand-alone test program.
187 * for example, compare
188 * a.out jupiter.9910 2451910.50000
189 * with
190 * satxyz2
191 * jup.dir.9910
192 * 2001 1 1 0 0 0
193 * 1 0 0 0
194 * 1
195 */
196int
197main (int ac, char *av[])
198{
199 double x[10], y[10], z[10];
200 char ynot[1024];
201 double jd;
202 char *fn;
203 FILE *fp;
204 int nm;
205 int i;
206
207 if (ac != 3) {
208 fprintf (stderr, "Usage: %s <filename> <jd>\n", av[0]);
209 abort();
210 }
211 fn = av[1];
212 jd = atof (av[2]);
213
214 fp = fopen (fn, "r");
215 if (!fp) {
216 perror (fn);
217 abort();
218 }
219
220 nm = read_bdl (fp, jd, x, y, z, ynot);
221 if (nm < 0) {
222 fprintf (stderr, "%s\n", ynot);
223 abort();
224 }
225
226 for (i = 0; i < nm; i++)
227 printf (" X= %19.11E Y= %19.11E Z= %19.11E\n", x[i], y[i], z[i]);
228
229 return (0);
230}
231#endif
232
233/* For RCS Only -- Do Not Edit */
234static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: bdl.c,v $ $Date: 2006-11-22 13:53:28 $ $Revision: 1.4 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.