Refaktorisierung von sweph()

Aufs Geratewohl herausgegriffen aus den verbliebenen knapp 5'000 Zeilen sei die Funktion sweph(). Ich stiess auf sie, weil sie von oben die erste Funktion ist, die die goto-Anweisung verwendet. Sicher ist goto considered harmful (Dijkstra) und sollte daher vermieden werden. Wichtiger ist aber die Verletzung des Single Responsibility Principle (SRP), dass eine Funktion nur für eine einzelne Aufgabe zuständig sein sollte. Das gilt ebenso für alle anderen Funktionen dieses Codes, aber irgendwo muss man einmal anfangen.

Die Funktion sweph() berechnet aus den Files der Swiss Ephemeris die baryzentrische Position des gewünschten Objekts. Gemäss ihrem einleitenden Kommentar hat sie folgende Aufgaben:

/*
 * this function looks for an ephemeris file,
 * opens it, if not yet open,
 * reads constants, if not yet read,
 * computes a planet, if not yet computed
 * attention: asteroids are heliocentric
 *            other planets barycentric
*/

Bei genauerer Inspektion des Codes führt sie die folgenden Aufgaben aus:

• Sie puffert den letzten Aufruf pro Planet. Die Struktur swed enthält für jeden Planeten einen "Slot" mit den Ergebnissen des letzten Aufrufs.
• Sie invalidiert den Zeiger auf das Ephemeris File, falls dieses von Datum oder Planetennummer nicht mehr gültig ist. In diesem Fall wird das File auch mit fclose() geschlossen und allokiertes Memory freigegeben.
• Ist das File ungültig, werden drei mögliche Filenamen für den korrekten Kandidaten ermittelt:
  1. Standardnamen wie sepl_18.se1 oder ast0/se000302.se1,
  2. sogenannte Kurzversionen mit einem "s" vor der File Extension wie in ast0/se00005s.se1,
  3. und schliesslich, falls der Pfad ein Verzeichnis enthielt wie in ast0/se000302.se1, wird ein letzter Versuch für den Filenamen ohne Verzeichnis gemacht, im Beispiel also se000302.se1.
  Die Kandidaten 2 und 3 werden dabei nur für Asteroiden probiert. Es wird versucht, Files mit diesem Namen zu öffnen. Das erste File, für das dies gelingt, wird das neue "aktuelle Ephemeris File".
• Es werden die für das ganze File gültigen Konstanten in die swed-Struktur eingelesen.
• Es wird am Anfang und Ende des Gültigkeitsbereichs der Swiss Ephemeris eine Randbedingung fürs Datum geprüft.
• Nun werden die Koeffizienten des für das Julianische Datum gültigen Interpolationspolynoms eingelesen.
• Das Polynom wird evaluiert.
• Für den Sonderfall der baryzentrischen Sonne wird die gesamte Funktion sweph() noch einmal für das Erde-Mond-Baryzentrum aufgerufen und das Ergebnis vektoriell zur Position der heliozentrischen Sonne addiert.
• Die Asteroidenpositionen sind heliozentrisch und werden in die baryzentrische Positionen konvertiert.
• Zum Kennzeichen, dass die berechneten Positionen gemerkt werden, wird die Zahl pdp->teval auf das Julianische Datum der Anfrage gesetzt.

Für eine Funktion mit 185 Zeilen sind das ziemlich viele Aufgaben. Auch gehören die Aufgaben zu verschiedenen Ebenen:

• Ebene File: Datei öffnen und schliessen, Filenamen
• Abbildung von Filedaten in Datenstruktur: Dateiweite Konstanten und Koeffizienten des zum aktuellen Objekt / Datum gültigen Interpolationspolynoms.
• Berechnung der baryzentrischen Position: Auswertung des Interpolationspolynoms, Umrechnung von Asteroidenpositionen von helio- in baryzentrisch.
• Pufferung der Abfrage

Nach Refaktorisierung dokumentiert der Code in der obersten Ebene die einzelnen Aufgaben, aus denen die Funktion besteht:

static int sweph(double tjd, int ipli, int ifno, int iflag, 
                 double *xsunb, AS_BOOL do_save, double *xpret, 
                 char *serr, struct swe_data *swed) {

  int ipl = is_main_object(ipli) ? ipli : SEI_ANYBODY;
  struct plan_data *pedp = &swed->pldat[SEI_EARTH];
  struct plan_data *psdp = &swed->pldat[SEI_SUNBARY];
  struct plan_data *pdp = &swed->pldat[ipl];

  /* In do_save mode, use swe_data memory slot for results */
  double* xp = do_save ? pdp->x : xpret;

  /* if planet has already been computed for this date, return */
  if (sweph_isPositionInBuffer(tjd, ipl, iflag, pdp, xpret))
    return (OK);

  /* Get ephemeris file handle */
  struct file_data* fdp = 
    sweph_getEphemerisFile(tjd, ipli, ifno, pdp, swed, serr);
  if (fdp == NULL )
    return NOT_AVAILABLE;

  if (sweph_isDateOutOfBound(tjd, fdp, serr, swed))
    return NOT_AVAILABLE;

  /* Read coefficient data */
  sweph_getCoefficientSegment(tjd, ipl, ifno, pdp, serr, swed);

  /* Also compute speed, if
   * 1. true position is being computed before applying light-time etc.
   *    this is the position saved in pdp->x.
   *    in this case, speed is needed for light-time correction.
   * 2. the speed flag has been specified.
   */
  bool need_speed = do_save || is_set(iflag, SEFLG_SPEED);

  /* Compute position */
  sweph_evaluateChebyshevPolynomial(tjd, pdp, need_speed, xp);

  /* Special case barycentric Sun */
  if (ipl == SEI_SUNBARY && (pdp->iflg & SEI_FLG_EMBHEL))
    sweph_barycentricSunFromEarth(tjd, ifno, iflag, need_speed, 
                                  pedp, xp, serr, swed);

  /* asteroids are heliocentric.
   * if JPL or SWISSEPH, convert to barycentric */
  if ((ipl >= SEI_ANYBODY) && 
      is_set(iflag, SEFLG_JPLEPH | SEFLG_SWIEPH))
    vec_add(xp, xsunb, need_speed ? 6 : 3);

  /* Mark result in save area as valid buffer */
  if (do_save) {
    sweph_validateBuffer(tjd, ifno, pdp, psdp);
    vec_copy(xpret, xp, 6);
  }

  return (OK);

}

Zwei Ebenen tiefer finden wir die zuvor mit goto again gelöste Aufgabe, das Ephemeridenfile anhand der drei Kandidaten für den Namen zu ermitteln. Die Code-Ersparnis, die das goto motiviert hatte, erwies sich als nicht besonders hoch: Man sparte sich nur die eine Codezeile mit dem Aufruf von swid_fopen(). Dadurch, dass es nun eine eigene kleine Funktkion, zwei Ebenen unterhalb von sweph() gibt, können wir return wie ein break benutzen, um auszusteigen, sobald swid_fopen() einen gültigen Filepointer lieferte:

static FILE* sweph_openEphemerisFile(int ipli, int ifno, char *fname, 
                                     char *serr, struct swe_data *swed) {
 
  FILE *fptr = swid_fopen(ifno, fname, swed->ephepath, serr, swed);
  if (fptr)
    return fptr;
 
  if (!is_main_object(ipli)) {
 
    char *sp;

    /*
     * if it is a numbered asteroid file:
     * try also for short files (..s.se1)
     */
    char s[AS_MAXCH];
    strcpy(s, fname);
    if ((sp = strrchr(s, '.'))) { 
      strcpy(sp, "s." SE_FILE_SUFFIX); /* insert an 's' */
      if ((fptr = swid_fopen(ifno, s, swed->ephepath, serr, swed)))
        return fptr;
      }
 
    /*
     * if we still have 'ast0' etc. in front of the filename,
     * we try in the main ephemeris directory instead of the
     * asteroid subdirectory.
     */
    if ((sp = strrchr(fname, *DIR_GLUE))) {
      if ((fptr = swid_fopen(ifno, sp+1, swed->ephepath, serr, swed)))
        return fptr;
    }
  
  }
  return NULL ;
}

Um festzustellen, wie Code für die Pufferung durchlaufen wird, schrieb ich einen neuen Test test_identical_call() für die Testsuite, der die Funktion se_calc zweimal nacheinander mit identischen Parametern aufruft. Überraschenderweise lieferte die Funktion sweph_isPositionInBuffer immer false!

Grund dafür war, dass exakt der gleiche Code, um festzustellen, ob das Ergebnis sich bereits im Puffer befindet, auch auf höherer Ebene bereits programmiert war - in der darüberliegenden Funktion sweplan(). Ich konnte diesen duplikativen Code nun durch Aufruf der Bedingungsfunktion ersetzen. Insgesamt gab es 6 identische Codestrecken in sweph.c, die jeweils diese Abfrage auf den Puffer machten, die ich nun alle durch den Aufruf der Bedingungsfunktion ersetzen konnte. Solche kleinen Bedingungsfunktionen werden übrigens vom Compiler meist automatisch ge-inlined. Die Bedenken, dass man durch den Funktionsaufruf eine zusätzliche Stackebene aufbauen und somit kostbare Rechenzeit verlieren würde, sind daher unbegründet.