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

Last change on this file since 3832 was 3654, checked in by cmv, 16 years ago

mise a niveau Xephem 3.7.4, cmv 16/07/2009
File size: 5.5 KB
RevLine 
[2551]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
[3654]117 if (width > sizeof(buf)-1) {
118 sprintf (ynot, "BDL Field width %d > %d", width, (int)sizeof(buf));
[2551]119 return (-1);
120 }
[3654]121 if (width != (int)fread (buf, 1, width, fp)) {
122 if (ferror(fp)) strcpy (ynot, "BDL IO error");
123 else if (feof(fp)) strcpy (ynot, "BDL unexpected EOF");
124 else strcpy (ynot, "BDL short file");
125 return (-1);
126 }
[2551]127
128 buf[width] = '\0';
129 switch (f) {
130 case I:
131 *(int *)ptr = atoi (buf);
132 break;
133 case F:
134 bp = strchr (buf, 'D');
135 if (bp)
136 *bp = 'e';
137 *(double *)ptr = atof (buf);
138 break;
139 case NL:
140 fgets (buf, sizeof(buf), fp);
141 break;
142 default:
143 sprintf (ynot, "Bug! format = %d", f);
144 return (-1);
145 }
146 return (0);
147}
148
149/* read one satellite record.
150 * return number of chars read else -1.
151 */
152static int
153readRec (FILE *fp, double *t0, double cmx[], double cfx[], double cmy[],
154double cfy[], double cmz[], double cfz[], char ynot[])
155{
156
157 long pos0, pos1;
158 int isat, idx;
159 int ldat1, ldat2;
160 int i;
161
162 pos0 = ftell (fp);
163
164 SCANFLD (I, 1, &isat);
165 SCANFLD (I, 5, &idx);
166 SCANFLD (I, 8, &ldat1);
167 SCANFLD (I, 8, &ldat2);
168 SCANFLD (F, 9, t0);
169 for (i = 0; i < 6; i++)
170 SCANFLD (F, 17, &cmx[i]);
171 for (i = 0; i < 4; i++)
172 SCANFLD (F, 17, &cfx[i]);
173 for (i = 0; i < 6; i++)
174 SCANFLD (F, 17, &cmy[i]);
175 for (i = 0; i < 4; i++)
176 SCANFLD (F, 17, &cfy[i]);
177 for (i = 0; i < 6; i++)
178 SCANFLD (F, 17, &cmz[i]);
179 for (i = 0; i < 4; i++)
180 SCANFLD (F, 17, &cfz[i]);
181 SCANFLD (NL, 0, NULL);
182
183 pos1 = ftell (fp);
184
185 return (pos1 - pos0);
186}
187
188
189#ifdef TEST_IT
190/* stand-alone test program.
191 * for example, compare
192 * a.out jupiter.9910 2451910.50000
193 * with
194 * satxyz2
195 * jup.dir.9910
196 * 2001 1 1 0 0 0
197 * 1 0 0 0
198 * 1
199 */
200int
201main (int ac, char *av[])
202{
203 double x[10], y[10], z[10];
204 char ynot[1024];
205 double jd;
206 char *fn;
207 FILE *fp;
208 int nm;
209 int i;
210
211 if (ac != 3) {
212 fprintf (stderr, "Usage: %s <filename> <jd>\n", av[0]);
213 abort();
214 }
215 fn = av[1];
216 jd = atof (av[2]);
217
218 fp = fopen (fn, "r");
219 if (!fp) {
220 perror (fn);
221 abort();
222 }
223
224 nm = read_bdl (fp, jd, x, y, z, ynot);
225 if (nm < 0) {
226 fprintf (stderr, "%s\n", ynot);
227 abort();
228 }
229
230 for (i = 0; i < nm; i++)
231 printf (" X= %19.11E Y= %19.11E Z= %19.11E\n", x[i], y[i], z[i]);
232
233 return (0);
234}
235#endif
236
237/* For RCS Only -- Do Not Edit */
[3654]238static char *rcsid[2] = {(char *)rcsid, "@(#) $RCSfile: bdl.c,v $ $Date: 2009-07-16 10:34:36 $ $Revision: 1.6 $ $Name: not supported by cvs2svn $"};
Note: See TracBrowser for help on using the repository browser.