Mittwoch, 21. Dezember 2011

Noch ein Zugang

Als Zwölfjähriger bekam ich zu Weihnachten ein Fahrrad geschenkt. Es war mir ein Vergnügen, dieses Fahrrad einmal vollständig bis zur letzten Schraube auseinanderzunehmen und danach wieder zusammenzubauen. Es blieben beunruhigenderweise am Schluss einige Schrauben übrig, aber der Fahrtüchtigkeit meines "Highrisers" tat das keinen Abbruch.

Diese Anlage - alles auseinanderzunehmen und wieder zusammenzubauen - könnte ich auch beim swepar-Projekt einbringen. Das ergibt - neben den vielen kleinen Refaktorisierungen, mit denen ich schon begonnen habe, einen weiteren Zugang zu diesem Projekt. Ich könnte damit beginnen, die Berechnung der baryzentrischen Koordinaten der Planeten aus der Swiss Ephemeris zu reproduzieren. Danach könnte das ganze Fleisch hinzugefügt werden: Umrechnung heliozentrisch -> geozentrisch, Berücksichtigung der Präzession, schliesslich Moshier als Ausweich-Ephemeride...

Die folgende kleine Funktion enthält alle nötigen Aufrufe, um die Position eines einzigen Planeten aus der Swiss Ephemeris zu berechnen. Das sieht doch eigentlich überschaubar aus...
 int test_direct() {

double tjd = 2452545.0;
int ipl = 5;
char qfname[255],fname[255],serr[255];
double xp[6];

// File name
swi_gen_filename(tjd, ipl, fname);
strcpy(qfname, "\\sweph\\ephe\\");
strcat(qfname, fname);

// Open the file
FILE* fp = fopen(qfname, BFILE_R_ACCESS);

// Put the minimum data necessary into swed
swed.fidat[0].fptr = fp;
strcpy(swed.fidat[0].fnam, fname);

// Read the necessary stuff
read_const(0, serr, &swed);

struct plan_data *pdp = & swed.pldat[ipl];
get_new_segment(tjd, ipl, 0, serr, &swed);
if (pdp->iflg & SEI_FLG_ROTATE)
rot_back(ipl,&swed);
else
pdp->neval = pdp->ncoe;

/* evaluate chebyshew polynomial for tjd */
double t = (tjd - pdp->tseg0) / pdp->dseg;
t = 2 * t - 1;
for (int i = 0; i <= 2; i++) {
xp[i] = swi_echeb (t, pdp->segp+(i*pdp->ncoe), pdp->neval);
xp[i+3] = swi_edcheb(t, pdp->segp+(i*pdp->ncoe), pdp->neval)
/ pdp->dseg * 2;
}

for (int i=0;i<6;i++) printf("xp[%d] = %15.10f\n",i,xp[i]);
return OK;
}