Klasse Mundan

package AstroOL;

/** Klasse für sphärisch-trigonometrische Berechnungen im Zusammenhang mit
dem <tt>motus primi mobilis</tt>: Umrechnungen zwischen Ekliptik und Äquator,
sowie von und in mundane Positionswerte. Zum Begriff der mundanen Position
vergleiche mein Buch "Primärdirektionen", Mössingen 1996, 4. Kapitel. */
public class Mundan {

/** Flag "Tag" (für Halbbogenberechnungen etc.) */
public static final byte day    = 0;
/** Flag "Nacht" (für Halbbogenberechnungen etc.) */
public static final byte night  = 1;
/** Flag "Tag und Nacht" (für Halbbogenberechnungen etc.) */
public static final byte always = 2;

/** Präzision bei Iterationen */
public static final double PREC = 1.e-8;
/** ungefaehr Bogensekundengenauigkeit */
public static final double  SPREC = 1.e-5;
/** Maximale Durchlaeufe bei Iterationen  */
public static final int  MAXIT = 10;
/** Maximale Zahl von Iterationen bei der Ber. der Eigenbewegung */
public static final int EBMAX = 3;

/** Häusersystem für mundane Positionsberechnungen */
public static final byte PLACIDUS      = 'P';
/** Häusersystem für mundane Positionsberechnungen */
public static final byte REGIOMONTANUS = 'R';
/** Häusersystem für mundane Positionsberechnungen */
public static final byte CAMPANUS      = 'C';
/** Häusersystem für mundane Positionsberechnungen */
public static final byte KOCH          = 'K';
/** Häusersystem für mundane Positionsberechnungen */
public static final byte PORPHYRIUS    = 'O';
/** Häusersystem für mundane Positionsberechnungen */
public static final byte AZIMUT        = 'H';
/** Häusersystem für mundane Positionsberechnungen */
public static final byte ALCABITIUS    = 'B';
/** Gesamtzahl implementierter Häusersysteme */
public static final int kMaxHouseSystem = 7;

/** Projektion längs der Breitenkreise (vertikal auf Ekliptik) */
public static final byte projectionVertical    = 0;
/** Projektion längs der Meridiane */
public static final byte projectionMeridional  = 1;
/** Projektion gemäss schräger Aufsteigung */
public static final byte projectionAscensional = 2;
/** Projektion längs mundaner Positionslinie */
public static final byte projectionMundan      = 3;


// Instanzvariablen
protected double lat = 0., st = 0., ee = 23.44;
protected byte hsys = PLACIDUS;

// Die häufigsten sin- und cos-Werte nur einmal ausrechnen
protected double sinEE = 0.397788507397949708158105598368606d,
                 cosEE = 0.917477140522918592843411380903052d,
                 tanLat = 0.,
                 sinSt = 0.;

/** Testmethode mit zufällig erzeugten Werten
    Läuft nur in VM 1.3 und höher
public static void main(String[] argv) {

  double a=0, d=0, l=0, b=0, mp=0, mp2;
  double[] x = new double[6];
  int i,iErr=0;
  byte system = PLACIDUS;
  boolean d0=false, l0=false, mp0=false;
  StringBuffer err = new StringBuffer();
  String houseSystems = "PRCKOHB";
  Mundan m;

  if (argv.length > 0) {

  for ( i = 3; i < argv.length; i++)
  switch (argv[i].charAt(0)) {
    case 'm':
       mp0 = true;
       mp = Double.parseDouble(argv[i].substring(1));
       break;
    case 'd':
       d0 = true;
       d = Double.parseDouble(argv[i].substring(1));
       break;
    case 'l':
       l0 = true;
       l = Double.parseDouble(argv[i].substring(1));
       break;
    case 's':
       system = Byte.parseByte(argv[i].substring(1));
    }

  m = new Mundan( Double.parseDouble(argv[0]),
                  Double.parseDouble(argv[1]),
                  Double.parseDouble(argv[2]),
                  system);



if (mp0) {
  if (d0) {
    a = m.radir(d,mp);
    m.cotrans(new double[] {a,d},x,true);
    l = x[0];
    b = x[1];
    mp = m.munpos(l,b,err);
    }
  else {
    l = m.hauslaenge(mp);
    b = 0;
    d = m.decl( l );
    a = m.ramc( l );
    mp = m.munpos(l,b,err);
    }
  }
else if (l0) {
  b = 0.;
  d = m.decl( l );
  a = m.ramc( l );
  mp = m.munpos(l,b,err);
  }

 System.out.println( "st\t"+Math.floor(1000*m.st)/1000d);
 System.out.println( "lat\t"+Math.floor(1000*m.lat)/1000d);
 System.out.println( "ee\t"+Math.floor(1000*m.ee)/1000d);
 System.out.println( "a\t"+Math.floor(1000*a)/1000d);
 System.out.println( "d\t"+Math.floor(1000*d)/1000d);
 System.out.println( "l\t"+Math.floor(1000*l)/1000d);
 System.out.println( "b\t"+Math.floor(1000*b)/1000d);
 System.out.println( "mp\t"+Math.floor(1000*mp)/1000d);
 System.out.println( "md\t"+Math.floor(1000*m.md(a, (mp>180)?day:night))/1000d);
 System.out.println( "hb\t"+Math.floor(1000*m.hb(d,(mp>180)))/1000d);
 System.out.println( "md/hb\t"+Math.floor(1000*m.md(a, (mp>180)?day:night)/m.hb(d,(mp>180)))/1000d);

 }

 else {

// Zufällige Tests: Breite zwischen -66 und +66,
// Sternzeit zwischen 0 und 360
// Deklination zwischen -40 und 40
// Mundanposition zwischen 0 und 360
// System 0...6

// radir(d,mp,system)     aufrufen   --> a
// munpos(a,d,system,serr) aufrufen  --> mp2

// arcdiffAbs(mp,mp2)>1.E-6 --> Error

  m = new Mundan();
  System.out.println("st\tlat\tsystem\tmp\tmp2\ta\td\tl\tb");
  for ( i = 0; i < 100; i++) {
    m.setSiderealTime( Math.random()*360 );
    m.setLatitude( Math.random()*132d - 66d );
    d = 80*Math.random() - 40;
    system = (byte)( Math.random() * kMaxHouseSystem);
    if (system == kMaxHouseSystem) system--;
    m.setHsys((byte) houseSystems.charAt(system));
    mp = Math.random()*360d;
    a = m.radir(d,mp);
    m.cotrans(new double[] {a,d},x,true);
    l = x[0];
    b = x[1];
    mp2 = m.munpos(l,b,err);
    if (err.length() > 0) System.out.println(err);
    if (Calculator.arcdiffAbs(mp,mp2)>1.E-6) {
       iErr++;
       System.out.print( Math.floor(1000*m.st)/1000d +"\t");
       System.out.print( Math.floor(1000*m.lat)/1000d +"\t");
       System.out.print( system +"\t");
       System.out.print( Math.floor(1000*mp)/1000d +"\t");
       System.out.print( Math.floor(1000*mp2)/1000d +"\t");
       System.out.print( Math.floor(1000*a)/1000d +"\t");
       System.out.print( Math.floor(1000*d)/1000d +"\t");
       System.out.print( Math.floor(1000*l)/1000d +"\t");
       System.out.println( Math.floor(1000*b)/1000d);
      }
    }
    System.out.println( iErr+" Fehler");



 }

}  */

/** Parameterloser Konstruktor. Vor Beginn der Rechnungen sind die
Sternzeit, geographische Breit und Ekliptikschiefe mit sinnvollen Werten
zu belegen! */
  public Mundan() {
  }

/** Explizite Konstruktion mit allen Parametern.
@param st Sternzeit (in Grad)
@param lat Geographische Breite
@param ee Ekliptikschiefe
@param hsys Häusersystem */
  public Mundan(double st, double lat, double ee, byte hsys ) {
    this.st = st;
    this.lat = lat;
    this.ee = ee;
    this.hsys = hsys;
    cosEE = Calculator.cos(ee);
    sinEE = Calculator.sin(ee);
    tanLat = Calculator.tan(lat);
    sinSt = Calculator.sin(st);
    }

/** Konstruktion durch ein Horoskop. Sternzeit, Breite und Ekliptikschiefe werden
dem Horoskop entnommen. */
  public Mundan( Horoscope h ) {
    this( h.st, h.lat, h.ee, h.hsys() );
    }

// public static Mundan getInstance(double xjd_ut, double lon, double lat, byte hsys) {
//  return null;
//  return new Mundan(sl.swe_sidtime( xjd_ut ) + lon, lat, sl.swi_epsiln(xjd_ut),hsys);
//  }

/** Setter-Methode für geographische Breite. */
public void setLatitude( double lat ) {
  this.lat = lat;
  tanLat = Calculator.tan( lat );
  }

/** Getter-Methode für geographische Breite. */
public double getLatitude() { return lat; }

/** Setter-Methode für Sternzeit (in Grad) */
public void setSiderealTime( double st ) {
  this.st = st;
  }

/** Getter-Methode für Sternzeit (in Grad) */
public double getSiderealTime() { return st; }

/** Setter-Methode für Häusersystem. */
public void setHsys( byte hsys ) {
  this.hsys = hsys;
  }

/** Getter-Methode für Häusersystem. */
public byte getHsys() {
  return hsys;
  }

/** Setter-Methode für Ekliptikschiefe. */
public void setEE( double ee ) {
  this.ee = ee;
  cosEE = Calculator.cos(ee);
  sinEE = Calculator.sin(ee);
  }

/** Getter-Methode für Ekliptikschiefe. */
public double ee() {
  return ee;
  }

/** Aszensionsdifferenz (AD, Unterschied von gerader und schräger Aufsteigung)
bei gegebener Deklination und Polhöhe.
@param d Deklination
@return Aszensionsdifferenz */
public static double ad(double d,double ph)
{
double x;
if (((x=Math.tan(d*Calculator.toRad)*Math.tan(ph*Calculator.toRad))<=1.)&&(x>=-1.))
 return Calculator.asin(x);
else {System.out.println("Fehler: Ascensionsdifferenz für Deklination " + d +
                         "und Breite " + ph + " undefiniert!"); return 0.;}
}

/** Aszensionsdifferenz (AD, Unterschied von gerader und schräger Aufsteigung)
bei gegebener Deklination und bei Polhöhe des Geburtsorts.
@param d Deklination
@return Aszensionsdifferenz */
public double ad( double d) { return ad( d, lat ); }

/** Deklination eines Ekliptikpunktes.
@param l Ekliptikale Länge
@return Deklination */
public double decl(double l) { return Calculator.asin(Calculator.sin(l)*sinEE); }

/** Deklination zu gegebener ekliptikaler Länge und Breite.
@param l Ekliptikale Länge
@param b Ekliptikale Breite
@return Deklination */
public double decl(double l, double b) {
  return Calculator.asin(
      Calculator.sin(l)*Calculator.cos(b)*sinEE +
      Calculator.sin(b)*cosEE);
  }

/** Schräge Aufsteigung zu gegebener Rektaszension, Deklination und Polhöhe.
@param a Gerade Aufsteigung (Rektaszension)
@param d Deklination
@param pol Polhöhe
@return Schräge Aufsteigung */
public static double oa(double a,double d, double pol)
{
return Calculator.fmod360(a-ad(d,pol));
}

/** Schräge Aufsteigung zu gegebener Rektaszension und Deklination (bei Polhöhe des Geburtsortes).
@param a Gerade Aufsteigung (Rektaszension)
@param d Deklination
@return Schräge Aufsteigung */
public double oa(double a,double d)
{
return Calculator.fmod360(a-ad(d,lat));
}

/** Schräge Aufsteigung zu gegebener ekliptikaler Länge (bei Polhöhe des Geburtsortes).
@param l Ekliptikale Länge
@return Schräge Aufsteigung */
public double oaasc(double l)
{
double y;
  y = ramc(l);
  y -= Math.asin(Math.sin(y*Calculator.toRad)*sinEE/cosEE*tanLat)*Calculator.toDeg; /* OAASC */
  return Calculator.fmod360(y);
}

/** Ekliptikale Länge zur Polhöhe des Geburtsortes und gegebener schräger Aufsteigung.
@param oa Schräge Aufsteigung
@return Ekliptikale Länge */
public double aclaenge(double oa) {  return aclaenge(oa,lat);  }


/** Ekliptikale Länge zu gegebener Polhöhe und schräger Aufsteigung.
@param oa Schräge Aufsteigung
@param phi Polhöhe
@return Ekliptikale Länge */
public double aclaenge(double oa, double phi) {
return Calculator.fmod360(Math.atan2(Math.sin(oa*Calculator.toRad),cosEE*Math.cos(oa*Calculator.toRad)-
   sinEE*Math.tan(phi*Calculator.toRad))/Calculator.toRad);
}

/** Rektaszension eines Ekliptikpunktes.
@param l Länge
@return Rektaszension */
public double ramc(double l)
{
return Calculator.fmod360(Math.atan2(Math.sin(l*Calculator.toRad)*cosEE,Math.cos(l*Calculator.toRad))*Calculator.toDeg);
}


/** Rektaszension zu gegebener ekliptikaler Länge und Breite.
@param l Länge
@param b Breite
@return Rektaszension */
public double ra(double l, double b ) {
  double[] coor = new double[2];
/* Umwandlung der ekliptikalen in aequatoriale Koordinaten */
  cotrans(new double[] {l,b},coor,false);
  return coor[0];
  }

/** Ekliptikale Länge zu vorgegebener Rektaszension.
@param a Rektaszension
@return Ekliptikale Länge */
public double mclaenge(double a)
{
return Calculator.fmod360(Math.atan2(Math.sin(a*Calculator.toRad),cosEE*Math.cos(a*Calculator.toRad))*Calculator.toDeg);
}

/** Meridionaldistanz. Distanz zum MC, reduziert auf das Intervall von 0 bis 360°.
@param a Rektaszension
@return Meridionaldistanz */
public double md( double a ) { return md(a, always); }

/** Meridionaldistanz zu MC oder IC.
@param a Rektaszension (in Grad)
@param visible <code>true</code> oder <code>false</code> für MC oder IC
@return Wert der Meridionaldistanz */
public double md( double a, boolean visible) { return md(a, visible? day:night); }

/** Meridionaldistanz zu MC oder IC.
@param a Rektaszension (in Grad)
@param tageszeit {@link #day day} oder {@link #night night} für MC oder IC; der
Wert {@link #always always} liefert die Distanz zum MC, auf das Intervall von 0 bis 360°
reduziert (im Gegensatz zu den beiden anderen Werten, bei denen der Wert auf das
Intervall (-180°;180°) reduziert wird, also um die Achse herum.
@return Wert der Meridionaldistanz */
public double md( double a, int tageszeit)
{
switch (tageszeit) {
  case day: return Calculator.arcdiff(a,st);
  case night: return Calculator.arcdiff(a,st+180.);
  case always: return Calculator.fmod360(a-st);
  default: System.out.println("Falscher Parameterwert (" + tageszeit + ") für Tageszeit!");
           return 0.;
  }
}

/** Halber Tag- oder Nachtbogen zur Deklination d.
@param d Deklination
@param tageszeit {@link #day day} oder {@link #night night} für Tag- oder Nachtbogen
@return Wert des Halbbogens */
public double hb( double d, byte tageszeit)
{
switch (tageszeit) {
  case always:
  case day: return 90.+ ad(d,lat);
  case night: return 90.-ad(d,lat);

  default: System.out.println("Falscher Parameterwert ("+tageszeit+") für Tageszeit!");
           return 0.;
           }
}

/** Halber Tag- oder Nachtbogen zur Deklination d.
@param d Deklination
@param visible Tag- oder Nachtbogen
@return Wert des Halbbogens */
public double hb( double d, boolean visible ) {
  return hb( d, visible?day:night );
  }

/** Steht der Himmelspunkt über oder unter dem Horizont? */
public boolean isVisible(double a, double d) {
  return Math.sin(ad(d,lat)*Calculator.toRad)+Math.cos((a-st)*Calculator.toRad)<0 ? false : true;
  }


/** Bestimmung der mundanen Position eines Himmelspunktes.
@param l     Ekliptikale Länge des Himmelspunktes
@param b     Ekliptikale Breite des Himmelspunktes
@return Mundane Position des Himmelspunktes.
*/
  public double munpos(double l, double b, StringBuffer serr)
{

int p;
double x,y,z,rc=0.,a,d;
double[] coor = new double[3];
boolean visible;

/* Voreinstellung: Keine Ambiguitaet */
if (serr != null) serr.setLength(0);

/* Umwandlung der ekliptikalen in aequatoriale Koordinaten */
cotrans( new double[] {l,b},coor,false);
a = coor[0];
d = coor[1];

switch (hsys) {

  case PLACIDUS:
    /* Prueft, ob Punkt ueber oder unter Horizont steht */

    visible = isVisible(a,d);
    if (visible) rc= 90.*(md(a,visible)/hb(d,visible)+3.);
    else rc= 90.*(md(a,visible)/hb(d,visible)+1.);
       break;

   case REGIOMONTANUS:
     if ((x=md(a,always))==0.) rc = 270.;
     /* Sonderfall obere Kulmination  */
     else if (x==180.) rc = 90.;
     /* Sonderfall untere Kulmination */
          else  rc = -Calculator.atan2(Calculator.tan(d)*tanLat+Calculator.cos(x),Calculator.sin(x));
          /* Allgemeine Formel */
     break;

   case CAMPANUS:
     if (d==90.||d==-90.) {System.out.println("Campanus ist für den Pol undefiniert!\n");
                           return 0.;}

     if ((x=md(a,always))==0) rc = 270.;
     /* Sonderfall obere Kulmination  */

     else if (x==180.) rc = 90.;
     /* Sonderfall untere Kulmination */
          else rc = -Calculator.atan2(
            Calculator.cos(lat)*Calculator.cos(d)*Calculator.cos(x)+
            Calculator.sin(lat)*Calculator.sin(d),Calculator.cos(d)*Calculator.sin(x));
          /* Allgemeine Formel */
     break;

   case KOCH:
      x=Calculator.asin(sinEE/cosEE*tanLat*sinSt);
      /* ADMC */
      y=ad(d); /* AD des Himmelspunktes */
      z=Calculator.arcdiff(a,st);
      /* MD als Wert zwischen -180 und +180 Grad */
      rc=0.;
      if (Calculator.arcdiffAbs(a+y,st+270.)<90.+x) /* WESTPUNKT */
        if (Calculator.arcdiffAbs(a+y,st+90)>=90.+x) /* & KEIN OSTPUNKT */
        {z+=180.;
         y=-y;
         rc=180.;
         }
         else
           if (b!=0) {
             serr.append("KOCH: Ambige Mundanposition. Punkt als Ostpunkt behandelt\n");
             System.out.println(serr);
             }
           else if ( Calculator.arcdiff(l,mclaenge(st)) < 0) {
             z+=180.;
             y=-y;
             rc=180.;
             }
      rc+=((z+x-y)/(90.+x)-1.)*90.;
    break;

   case PORPHYRIUS:
     x = aclaenge(st+90.);
     y = mclaenge(st);
     if ( (a = Calculator.arcdiff(x,l))  >= 0)
       if ( (z = Calculator.arcdiff(y,l)) >= 0)
         rc =  270. - 90.*z/Calculator.arcdiff(y,x-180.);  /* Dritter Quadrant */
       else
         rc =  270. - 90.*z/Calculator.arcdiff(x,y);       /* Vierter Quadrant */
     else
       if ( (z = Calculator.arcdiff(y+180.,l)) >= 0)
         rc = 90. - 90.*z/Calculator.arcdiff(y+180.,x);    /*  Erster Quadrant */
       else
         rc = 90. - 90.*z/Calculator.arcdiff(x,y);         /* Zweiter Quadrant */
     break;

   case AZIMUT:
     p = Calculator.sgn( tanLat-sinEE/cosEE*sinSt );
     if (lat == 0) return Calculator.atan2(p*Calculator.sin(d),Calculator.cos(d)*Calculator.sin(a-st));
     if ( (x = Calculator.sin(a-st)) == 0)
       if (a == st) return (p*(lat-d)>=0) ? 270.:90.;
       else return (p*(lat-d)<0) ? 270. : 90.;
     else {
       y = Calculator.tan(d)/Calculator.fabs(tanLat)-Calculator.sgn(lat)*Calculator.cos(a-st);
       rc = p*Math.atan((Calculator.cos(lat)*Calculator.tan(d)
                         -Calculator.cos(a-st)*Calculator.sin(lat))/x)*Calculator.toDeg;
       if (rc*p*y<0) rc +=180.;
       }
     break;

   case ALCABITIUS:
     x = Calculator.asin(Calculator.sin(aclaenge(st+90.))*sinEE);
     /* Quadranten ermitteln */
     p = (Calculator.arcdiff(a,st + hb(x,day)) > 0 ) ? 0 : 2; /* Nachtpunkt? */
     p += (Calculator.arcdiff(a,st + 90.*(p-2)) > 0) ? 2 : 1;         /* Westpunkt?  */
     switch (p) {
       case 1:
         rc = (1. + Calculator.arcdiff(a,st+180.)/hb(x,night))*90.;
         break;

       case 2:
         rc = (1. + Calculator.arcdiff(a,st+180.)/hb(x,day))*90.;
         break;

       case 3:
         rc = (3. + Calculator.arcdiff(a,st)/hb(x,night))*90.;
         break;

       case 4:
         rc = (3. + Calculator.arcdiff(a,st)/hb(x,day))*90.;
     }
     break;

   default: serr.append("Ambige Mundanposition. Punkt als Ostpunkt behandelt\n");
            System.out.println(serr);
            return 0.;
  }  /* switch */
return Calculator.fmod360(rc);
} /* munpos */


/** Ekliptikale Länge zu gegebener mundaner Position.
@param mp Mundane Position
@return Länge des Ekliptikpunktes der mundanen Position mp. Für mp = 0, 30,
60, usw. sind dies die Längen der Häuserspitzen. */
public double hauslaenge( double mp)
{
double alpha,a,c,x=0.,y,rc=0.;
short p;

switch (hsys) {

case PLACIDUS:
  alpha  = Calculator.fmod360(st+mp+90.)*Calculator.toRad;
  /* In Bogenmass umgerechnete OA als Startwert */
  x=alpha; /* OA */
  a= (Calculator.arcdiff(mp,180)>0) ? (mp-270.)/90. : (90.-mp)/90.;
  c = sinEE/cosEE*tanLat;

  /* Iteration: */
  do  alpha -= (y=(x+a*Math.asin(c*Math.sin(alpha))-alpha) / (a*c*Math.cos(alpha)-1));
    while (Calculator.fabs(y)>PREC);

  return mclaenge(alpha/Calculator.toRad);
  /* Ergebnis-RA in ekliptikale Laenge verwandeln */

case REGIOMONTANUS:
  x=(st+mp);
  y=mp;
  return Calculator.fmod360(-Calculator.atan2(Calculator.cos(x),cosEE*Calculator.sin(x)+sinEE*tanLat*Calculator.cos(y))+180.);

case CAMPANUS:
  x=Calculator.asin(Calculator.cos(mp)*Calculator.sin(lat));                             /* Polhoehe */
  y=-Calculator.atan2(Calculator.cos(mp)*Calculator.cos(lat),Calculator.sin(mp))+st+180.;                  /* OA */
  return aclaenge(y,x);
                 /* Ekliptikale Laenge zu dieser OA und Polhoehe */

case KOCH:
   y=x=Calculator.fmod360(mp+90.)/90.-1.; /* auf -90...+270 reduzierte mp */
   c=lat;  /* Polhoehe */
   if (x>1) {x-=2;c=-lat;}   /* Bei Westpunkten: Mit OD arbeiten */
   y=st+90.*(1+y) + x*Calculator.asin(sinEE/cosEE*tanLat*sinSt);           /* OA */
   return aclaenge(y,c);
  /* Ekliptikale Laenge zu dieser OA und Polhoehe                */

case PORPHYRIUS:
  x = aclaenge(st+90.);
  y = mclaenge(st);
  switch ((short) (mp/90.)) {
    case 0:
      rc = x + Calculator.arcdiff(y+180.,x)*mp/90;
      break;
    case 1:
      rc = y + 180. + Calculator.arcdiff(x,y)*(mp/90.-1);
      break;
    case 2:
      rc = x + 180. + Calculator.arcdiff(y+180.,x)*(mp/90.-2);
      break;
    case 3:
      rc = y + Calculator.arcdiff(x,y)*(mp/90.-3);
    }
  return Calculator.fmod360(rc);

case AZIMUT:
  p = Calculator.sgn(tanLat-sinEE/cosEE*sinSt);
  x = -Calculator.sgn(lat)*Calculator.asin(Calculator.cos(mp)*Calculator.cos(lat));
  y = -Calculator.asin(tanLat*Calculator.tan(x));
  if (mp<180.) y = 180. - y;
  return aclaenge(y+st,x);

case ALCABITIUS:
  c = Calculator.asin(sinEE*Calculator.sin(aclaenge(st+90.)));
  y = Calculator.fmod360(mp)/90.;
  switch ( (short) y + 1) {   /* Quadranten */
    case 1:
      x = (y-1.)*hb(c,night) + st + 180.;
      break;

    case 2:
      x = (y-1.)*hb(c,day) + st + 180.;
      break;

    case 3:
      x = (y-3.)*hb(c,night) + st;
      break;

    case 4:
      x = (y-3.)*hb(c,day) + st;
    }
  rc = mclaenge(x);
  return rc;

default: {System.out.println("Hauslaenge: Haeusersystem" +hsys+ "unbekannt!"); return 0;}

}  /* SWITCH */

}

/** Ekliptikale Länge zu gegebener Sternzeit und mundaner Position.
@param mp Mundane Position
@param st Sternzeit
@return Länge des Ekliptikpunktes der mundanen Position mp. Für mp = 0, 30,
60, usw. sind dies die Längen der Häuserspitzen. */
public double hauslaenge(double mp,double st) {
  double st_save = this.st,l;
  setSiderealTime( st );
  l = hauslaenge( mp );
  setSiderealTime( st_save );
  return l;
  }

/** Rektaszension zu gegebener Deklination und mundaner Position.
@param d      Deklination
@param mp     Mundane Position
@return Rektaszension des Punktes auf dem durch d gegebenen
Parallelkreis, der die mundane Position mp hat. */
public double radir(double d,double mp)

{
double a,x=0,y,c;
short p;

switch (hsys) {

case PLACIDUS:
  if (Calculator.arcdiff(mp,180.)>=0) a=(mp-270.)/90.;
  else a=(90.-mp)/90.;
  return Calculator.fmod360(mp+st+a*ad(d)-270.);
  /* Zu OA den Anteil der Aszensionsdifferenz aufaddieren */

case REGIOMONTANUS:
  return Calculator.fmod360(Calculator.acos(-Calculator.tan(d)*tanLat*Calculator.cos(mp))+st+mp);

case CAMPANUS:
  x=Calculator.asin(Calculator.cos(mp)*Calculator.sin(lat));            /* Polhoehe */
  y=-Calculator.atan2(Calculator.cos(mp)*Calculator.cos(lat),Calculator.sin(mp))+st+180.; /* OA */

  return Calculator.fmod360(y+ad(d,x));  /* RA = OA + AD */

case KOCH:
   x=Calculator.fmod360(mp+90.)/90.-1.;   /* auf -90...+270 reduzierte mp */
   c=lat;  /* Polhoehe */

   if (x>1) {x-=2;c=-lat;}  /* Bei Nachtpunkten: Mit OD arbeiten */

   y=st+90.*(1+x) + x*Calculator.asin(sinEE/cosEE*tanLat*sinSt);  /*          OA */

   return Calculator.fmod360(y+ad(d,c));                  /* RA = OA + AD */

case PORPHYRIUS:
  x = hauslaenge(mp);
  c = Calculator.atan2(Calculator.sin(x),Calculator.cos(x)*cosEE);
  y = Calculator.asin(sinEE*Calculator.cos(x));
  return Calculator.fmod360( c-Calculator.asin(Calculator.tan(d)*Calculator.tan(y)));

case AZIMUT:
  x = Calculator.asin( -Calculator.sgn(lat)*Calculator.cos(mp)*Calculator.cos(lat));
  y =  - Calculator.asin(Calculator.tan(x)*tanLat);
  if (mp <180.) y = 180. - y;
  return Calculator.fmod360(y+st+Calculator.asin(Calculator.tan(d)*Calculator.tan(x)));

case ALCABITIUS:
  c = Calculator.asin(sinEE*Calculator.sin(aclaenge(st+90.)));
  y = Calculator.fmod360(mp)/90.;
  switch ( (short) y + 1) {   /* Quadranten */
    case 1:
      x = (y-1.)*hb(c,night) + st + 180.;
      break;

    case 2:
      x = (y-1.)*hb(c,day) + st + 180.;
      break;

    case 3:
      x = (y-3.)*hb(c,night) + st;
      break;

    case 4:
      x = (y-3.)*hb(c,day) + st;
    }
  return Calculator.fmod360(x);

default:
   System.out.println("Radir: Unbekanntes Haeusersystem: "+hsys);
   return 0.;

} /* SWITCH */
} /* RADIR */

  /** Die Routine <code>swe_cotrans()</code> der Swiss Ephemeris
  kann aus technischen Gründen nicht eingesetzt werden, da die
  Klasse Mundan auch in Applets verwendet wird und das
  volle Swisseph-Paket zu grosse Ladezeiten hervorrufen würde.
   */
  public void cotrans(double in[], double out[], boolean toEcl) {
    out[0] = Calculator.atan2(
      Calculator.sin(in[0])*cosEE + Calculator.tan(in[1])*sinEE*(toEcl?1:-1),
      Calculator.cos(in[0]));
    if (out[0] < 0) out[0] += 360d;
    out[1] = Calculator.asin(
      Calculator.sin(in[0])*Calculator.cos(in[1])*sinEE*(toEcl?-1:1) +
      Calculator.sin(in[1])*cosEE);

  }

}