Compute a transit list

// Compute a list of transits over aspects of a given horoscope
// Also compute the enter and leave times into the 1° orb of each aspect

#include <iostream>
#include <iomanip>
#include <vector>
#include <regex>
#include <map>
#include <math.h>

#include "swephexp.h"

struct t_input {
  double jd_et,lon,lat;
  double from,to;
  unsigned int filter;
};

using namespace std;

// Difference on the arc (taking account of the gap 360->0)
inline double arcdiff(double x, double y) {
  double z = x - y;
  return ( (z >= -180) && (z < 180)) ? z :
    z - floor(z/360+.5)*360;
}

// Reduce a degree value to the [0...360[
inline double fmod360( double x ) {
  return (x < 360 && x >= 0) ? x :
    x - floor(x/360)*360;
}

// An ecliptical position (only longitude may be filled)
// Also keeps track of possible warnings during computation
class Position {
  public:
    Position() : Position(0.) {}
    Position(const double lon, const string& warning = "") {
      this->lon   = lon;
      this->lat   = 0;
      this->dist  = 0;
      this->vlon  = 0;
      this->vlat  = 0;
      this->vdist = 0;
      this->warning = warning;
    }
    Position(const double *xx, const string& warning ="") {
      this->lon   = xx[0];
      this->lat   = xx[1];
      this->dist  = xx[2];
      this->vlon  = xx[3];
      this->vlat  = xx[4];
      this->vdist = xx[5];
      this->warning = warning;
    }
    double lon, lat, dist;
    double vlon,vlat,vdist;
    string warning;
    bool has_warnings() {
      return !warning.empty();
    }
    void toStream(ostream &str) const {
      str << lon;
    }
};

// Ouput of a position
ostream & operator<<(ostream & Str, Position const & p) {
  p.toStream(Str);
  return Str;
}

class ComputationError {
public:
  ComputationError(const string message) : message(message) {}
  string message;
};

// The minimum which is necessary for computing a horoscope: jd, lon, lat
class HoroscopeBasicData {
public:
  double jd_et,lon,lat;
  bool initial = true;
  HoroscopeBasicData() {}
  HoroscopeBasicData(double jd_et,double lon, double lat) : jd_et(jd_et),lon(lon),lat(lat) {
    initial = false;
  }
  bool operator == (const HoroscopeBasicData& h) {
    return h.jd_et == jd_et && h.lon == lon && h.lat == lat;
  }
  bool operator != (const HoroscopeBasicData& h) {
    return h.jd_et != jd_et || h.lon != lon || h.lat != lat;
  }
};


// An astrological object
// abstraction of all computable points of a Horoscope
// (like planets, asteroids, fixed stars, house cusps, midpoints, ...)
class AstroObject {

public:
  const string name;
  const string symbol;
  unsigned int aspectStrength = 0x7f; // Default: aspects of all groups

  AstroObject(const string& name, const string& symbol) :
     name(name),symbol(symbol) {}
  AstroObject(const char* name, const char* symbol) : name(name), symbol(symbol) {}

  virtual Position compute( const HoroscopeBasicData& data ) const throw(ComputationError) = 0;
  virtual ~AstroObject() {}

protected:
  friend ostream& operator<<(ostream&, const AstroObject&);

};
ostream& operator<<(ostream &strm, const AstroObject &a) {
  return strm << a.symbol ;
}

// I define an ephemeris object by computability via a swe_calc() call
class EphemerisObject : public AstroObject {
public:
  EphemerisObject(const string& name, const string& symbol, const int number):
    AstroObject(name,symbol), ipl(number) {
  }
  Position compute( const HoroscopeBasicData& data ) const throw(ComputationError) {
    double xx[6];
    char serr[256];
    int iflag = SEFLG_SPEED;
    int rc = swe_calc( data.jd_et, ipl, iflag, xx, serr);
    if (rc < 0) {
      throw ComputationError(string(serr));
    }
    return Position(xx,rc > 0 ? string(serr) : "");
  }
private:
  int ipl;
};


class Planet : public EphemerisObject {
public:
  Planet(const string& name, const string& symbol, const int number):
    EphemerisObject(name,symbol,number) {
  }
};

// We want to treat a single house cusp like any other AstroObject
// On the other hand, swe_houses gives us all the house cusps at once
// For efficiency, we have a class Houses treating the complete group...
class Houses;

// ... and a class House representing a single house cusp
class House : public AstroObject {
public:
  House(const char* name, const char* symbol, const int number, Houses *houses):
    AstroObject(name,symbol), number(number), houses(houses) {
    // Intermediate house cusps: ☌☌ only
    // Angles: ☌☌, □□, △△ (respect symmetry of AC/DC and MC/IC)
    switch (number) {
      case 1:
      case 10:
       aspectStrength = 0b1101;
       break;
      default:
       aspectStrength = 0b1;
    }
  }
  Position compute( const HoroscopeBasicData& data ) const throw(ComputationError);
private:
  int number;
// Reference to a Houses instance, which computes all the house cusps
  Houses* houses;

};


class Houses {
public:
  struct t_cusp {
    const char * name, *symbol;
  };
  Houses(int hsys = 'P',  // Default: Placidus
      const t_cusp cusps[] = (t_cusp[]) {
      {"ASC","AC"},
      {"Haus 2","H2"},
      {"Haus 3","H3"},
      {"IC","IC"},
      {"Haus 5","H5"},
      {"Haus 6","H6"},
      {"DSC","DC"},
      {"Haus 8","H8"},
      {"Haus 9","H9"},
      {"MC","MC"},
      {"Haus 11","H11"},
      {"Haus 12","H12"},
      {NULL,NULL}
    }) : hsys(hsys) {
    for (int i=0;cusps[i].name != NULL;i++) {
      this->cusps.push_back( new House( cusps[i].name, cusps[i].symbol, i+1, this) );
    }
  }
  vector<const AstroObject*> cusps;
  int hsys;

  map<const AstroObject*,Position> compute( const HoroscopeBasicData& data ) throw(ComputationError);
  Position compute( const HoroscopeBasicData& data, const int number ) throw(ComputationError);

private:
  HoroscopeBasicData data;
  vector<Position> positions;

};

// Default values: 12 house cusps, Placidus
Houses* housesDefault = new Houses();

// Default: the ten most commonly used astrological planets, plus lunar node
vector<const AstroObject*> planetsDefault = {
  new Planet("Sonne","☉",0),
  new Planet("Mond","☽",1),
  new Planet("Merkur","☿",2),
  new Planet("Venus","♀",3),
  new Planet("Mars","♂",4),
  new Planet("Jupiter","♃",5),
  new Planet("Saturn","♄",6),
  new Planet("Uranus","⛢",7),
  new Planet("Neptun","♆",8),
  new Planet("Pluto","♇",9),
  new Planet("Mondknoten","☊",10)
};


// Descriptive class, only keeps a vector of AstroObjects, in particular house cusps
// Has a compute() method for getting all the positions in a concrete Horoscope
class HoroscopeParts {
public:
  HoroscopeParts( vector<const AstroObject*>& planets, Houses* houses ) : planets(planets),houses(houses) {}
  map<const AstroObject*,Position> compute( const HoroscopeBasicData& data ) const {
    map<const AstroObject*,Position> result;
    for (const AstroObject* p : planets ) {
      result[p] = p->compute( data );
    }
    map<const AstroObject*,Position> hpos = houses->compute( data );
    result.insert( hpos.begin(), hpos.end() );
    return result;
  }
private:
  vector<const AstroObject*> planets;
  Houses* houses;
};

HoroscopeParts horoscopePartsDefault( planetsDefault, housesDefault);

// Implements the swe_houses() call
map<const AstroObject*,Position> Houses::compute(const HoroscopeBasicData& data) throw (ComputationError) {
  map<const AstroObject*,Position> result;
  this->data = data;
  positions.clear();
  double houses[13], ascmc[10];
  swe_houses( data.jd_et - swe_deltat( data.jd_et), data.lat, data.lon, hsys, houses, ascmc);
  for (int i=0;i<12;i++) {
    Position p = Position( houses[i+1] );
    result[ cusps[i] ] = p;
    positions.push_back( p );
  }
  return result;
}

// Give back a single house cusp position
Position Houses::compute( const HoroscopeBasicData& data, const int number ) throw(ComputationError) {
    if (this->data.initial ||
        data.jd_et != this->data.jd_et ||
        data.lon   != this->data.lon   ||
        data.lat   != this->data.lat ) {
      compute( data );
    }
    return positions[number-1];
}

// The same, viewed from inside of AstroObject
Position House::compute( const HoroscopeBasicData& data ) const throw(ComputationError) {
  return houses->compute(data,number);
}

// Descriptive class, objects representing single aspects
class Aspect {
public:
  Aspect(
    const char* name,
    const char* symbol,
    int strength,
    double angle) : name(name), symbol(symbol), strength(strength), angle(angle) {}
  string name;
  string symbol;
  unsigned int strength;
  double angle;
  ostream& toStream( ostream& str ) const {
    str << symbol;
    return str;
  }
};

ostream & operator<<(ostream & str, Aspect const & a) {
  return a.toStream(str);
}

// The aspects in use, and their groups
vector<Aspect> aspects = {
  { "Konjunktion",        "☌",  0x01, 0  },
  { "Halbsextil",         "⚺",  0x40, 30 },
  { "Halbquadrat",        "∠",  0x20, 45 },
  { "Sextil",             "⚹",  0x10, 60 },
  { "Quadrat",            "□",  0x04, 90 },
  { "Trigon",             "△",  0x08, 120 },
  { "Anderthalbquadrat",  "⚼",  0x20, 135},
  { "Quincunx",           "⚻",  0x40, 150},
  { "Opposition",         "☍",  0x02, 180}
};

enum Direction {
  LEFT,
  RIGHT
};

// Descriptive structure for an aspect to a planet (to an AstroObject)
class AspectPoint {
public:
  const Aspect* aspect;       // Aspect...
  const AstroObject* object;  // ...to this AstroObject
  const Direction direction;  // ...RIGHT = in forward, LEFT = in backward zodiac direction
  const int bound;            // Usually zero, but +1 or -1 for aspect boundaries

};

ostream & operator<<(ostream & str, AspectPoint const & ap) {
  str << *(ap.aspect) << *(ap.object);
  return str;
}


// Compputes an aspect wreath from a map of AstroObjects with positions
class Aspectarium {
public:
  multimap<double,AspectPoint> wreath;
  Aspectarium(const map<const AstroObject*,Position> p, const unsigned int filter = 0xff, const vector<Aspect>& aspects = aspects) {
    for (auto o : p) {
      for (const Aspect& a: aspects) {
        if (!(a.strength & filter )) continue;
        if (!(o.first->aspectStrength & a.strength)) continue;
        addAspect(o.first,o.second.lon,RIGHT,a);
        if (a.symbol != "☌" && a.symbol != "☍") {
          addAspect(o.first,o.second.lon,LEFT,a);
        }
      }
    }
  }

  multimap<double,AspectPoint> wreath_enriched_with_boundaries() const {
    multimap<double,AspectPoint> result = wreath;
    for (auto it : wreath) {
      result.insert( pair<double,AspectPoint>(
        fmod360(it.first-1),
        AspectPoint({ it.second.aspect,it.second.object,it.second.direction,-1 })
      ));
      result.insert( pair<double,AspectPoint>(
        fmod360(it.first+1),
        AspectPoint({ it.second.aspect,it.second.object,it.second.direction,+1 })
      ));
    }
  return result;
  }

  ostream& toStream( ostream& str) const {
    for (auto a: wreath) {
      AspectPoint ap = a.second;
      str << setw(12) << setprecision(4) << fixed << a.first << " : " << ap << endl;
    }
    return str;
  }
private:
  void addAspect(const AstroObject* o, const double pos, const Direction direction, const Aspect& a) {
    double l = fmod360( pos + a.angle * (direction == RIGHT ? 1 : -1) );
    AspectPoint ap({ &a, o, direction, 0 });
    wreath.insert( pair<double,AspectPoint>(l,ap));
  }

};

ostream & operator<<(ostream & str, Aspectarium const & a) {
  a.toStream(str);
  return str;
}


// Class Horoscope, basically wraps a map of AstroObjects together with their positions
class Horoscope : public HoroscopeBasicData {
public:
  Horoscope(double jd_et,double lon, double lat,const HoroscopeParts& hp = horoscopePartsDefault) : HoroscopeBasicData(jd_et,lon,lat) {
    p = hp.compute((HoroscopeBasicData)*this);
  }

  Aspectarium getAspectarium( unsigned int aspect_filter = 0xff) {
    return Aspectarium( p, aspect_filter );
  }

  void toStream(ostream &str) const {
    for (auto const& pos: p) {
      str << setw(20) << *(pos.first) << ": " << pos.second << endl;
    }
  }

private:
  map<const AstroObject*,Position> p;

};

ostream & operator<<(ostream & Str, Horoscope const & h) {
  h.toStream(Str);
  return Str;
}

// Position with time
typedef struct {
  double jd_et;
  Position pos;
} t_position;

// Any interval
typedef struct {
  double from, to;
} t_interval;

// Compute the time of a transit
double transit_time(
  double lon,
  const t_position last,
  const t_position next) {

  double jd_et;

  if (next.pos.lon == last.pos.lon) return (last.pos.lon == lon) ? last.jd_et : NAN;

// To keep it simple, we use linear interpolation
// Could improve by cubic interpolation, employing also last.pos.vlon and next.pos.vlon
  jd_et =
    (next.jd_et-last.jd_et)*arcdiff(lon,last.pos.lon)/arcdiff(next.pos.lon,last.pos.lon)
    + last.jd_et;

  return jd_et;

}

// An object of class Transit represents a single transit
// a current AstroObject is crossing an aspect to a given AstroObject (from the base horoscope)
class Transit {
public:
  double jd_et,lon;
  const AstroObject* current;
  AspectPoint aspect;
  Direction direction;
  double jd_begin, jd_end;


  // Current planet is entering the orbis
  // - either it is moving directly, and the bound is -1
  // - or it is moving retrogradely, and the bound is +1
  bool entering() const {
    return (
      (aspect.bound == -1 && direction == RIGHT) ||
      (aspect.bound == +1  && direction == LEFT)
    );
  }

  // Current planet is leaving the orbis
  // - either it is moving directly, and the bound is +1
  // - or it is moving retrogradely, and the bound is -1
  bool leaving() const {
    return (
      (aspect.bound == -1 && direction == LEFT) ||
      (aspect.bound == +1  && direction == RIGHT)
    );
  }

  // Consider two transit objects as equal if the aspect formulae are identical
  // here, we do distinguish right from left aspects
  bool operator==(const Transit& t) const {
    return (
     current       == t.current &&
     aspect.aspect == t.aspect.aspect &&
     aspect.object == t.aspect.object &&
     aspect.direction == t.aspect.direction );
  }
};

// Output a transit list as JSON object
void print(ostream & Str, const multimap<double,Transit>& t) {
  bool first = true;
  Str << "[" << endl;
  for ( auto it : t ) {
    if (!first) {
      Str << "," << endl;
    }
    Str << "["
         << setw(14) << setprecision(5) << fixed << it.first << ", "
         << setw(14) << setprecision(5) << fixed << it.second.jd_begin << ", "
         << setw(14) << setprecision(5) << fixed << it.second.jd_end << ", "
         << setw(8) << setprecision(4) << it.second.lon << ", "
         << "\"" << *(it.second.current)
         << (it.second.direction == LEFT ? "℞" :"")
         << (it.second.aspect.bound == -1 ? "LB " :"")
         << (it.second.aspect.bound == 1 ? "RB " :"")
         << it.second.aspect << "\"]";
    first = false;
  }
  Str << endl << "]";
}

// Output an aspect wreath as a list (plain text, not JSON)
void print(ostream & Str, const multimap<double,AspectPoint>& w) {
  for (auto it: w) {
    Str << setw(14) << setprecision(5) << fixed << it.first << ": ";
    Str << it.second << endl;
  }
}

// Find all the aspects from the wreath falling into the given interval,
// and compute the exact transit time for each of them
multimap<double,Transit> transits_from_interval(
  const t_interval& interval,
  const AstroObject* o,
  const t_position last,
  const t_position next,
  const multimap<double,AspectPoint>& wreath) {

  multimap<double,Transit> r;
  for (auto it = wreath.lower_bound( interval.from );
       it != wreath.end() && it->first < interval.to;
       ++it ) {
    double jd_et = transit_time(it->first,last,next);
    r.insert( pair<double,Transit>(
      jd_et,
      Transit({
        jd_et,
        it->first,
        o,
        it->second,
        (arcdiff(next.pos.lon,last.pos.lon)<0) ? LEFT: RIGHT })
    ));
  }
  return r;
}

// Decompose an interval on the circle into parts on [0...360]
// Take account of the discontinuity 360->0
// Take account of retrograde motion (input.from may well be > input.to)
vector<t_interval> get_segments(const t_interval& input) {
  if (arcdiff(input.to,input.from)>0) {
    if (input.to > input.from) {
      return vector<t_interval>({input});
    } else {
      return vector<t_interval>({
        { input.from, 360.},
        { 0., input.to }
      });
    }
  } else {
    if (input.to < input.from) {
      return vector<t_interval>({{input.to,input.from}});
    } else {
      return vector<t_interval>({
        { input.to, 360.},
        { 0., input.from}
      });
    }
  }
}

// Compute all transits in the given interval,
// including the times when the aspect is entered and left by the current AstroObject
multimap<double,Transit> transits(double from, double to, const Aspectarium& a, vector<const AstroObject*>& planets = planetsDefault) {

// Enrich the aspect wreath with aspect boundaries (i.e. aspect locus +1° / -1°)
  multimap<double,AspectPoint> wreath = a.wreath_enriched_with_boundaries( );

// For each planet:
// Compute the times where it passes the aspects of the 'enriched wreath'
  multimap<double,Transit> r;
  for (const AstroObject* o : planets) {
    Position last = o->compute( HoroscopeBasicData(from,0,0) );
    for (double t = from+1; t<=to; t+=1) {
      Position next = o->compute( HoroscopeBasicData(t,0,0) );
      vector<t_interval> v = get_segments({ last.lon, next.lon });
      for (t_interval iv: v) {
        multimap<double,Transit> s = transits_from_interval(
          iv,
          o,
          (t_position) { t-1, last },
          (t_position) { t, next },
          wreath
        );
        r.insert(s.begin(),s.end());
      }
      last = next;
    }
  }

// For each transit: adjoin the time when it enters and leaves the orbis
  multimap<double,Transit> result;
  for (auto it = r.begin(); it != r.end(); ++it) {
    if (it->second.aspect.bound == 0) {

// Scan backwards in longitude, when did the planet enter the aspect orb
      it->second.jd_begin = from;
      for (auto it1 = it; it1 != r.begin();--it1) {
        if (it1->second == it->second &&
            it1->second.entering()) {
          it->second.jd_begin = it1->first;
          break;
        }
      }

// Scan forward in longitude, when did the planet leave the aspect orb
      it->second.jd_end = to;
      for (auto it2 = next(it); it2 != r.end();++it2) {
        if (it2->second == it->second &&
            it2->second.leaving() ) {
          it->second.jd_end = it2->first;
          break;
        }
      }

      result.insert(*it);
    }
  }
  return result;
}



void teardown() {

  for (auto p: planetsDefault) {
    delete p;
  }

}

void parse_query_string( const char* qs, t_input& input ) {

    const regex NAME_VALUE_PATTERN("([\\w+%]+)=([^&]*)");
    const string query(qs);
    auto words_begin = sregex_iterator(query.begin(), query.end(), NAME_VALUE_PATTERN);
    auto words_end = sregex_iterator();

    for (auto it = words_begin; it != words_end; it++) {
      const string name((*it)[1]);
      const char* value = (*it)[2].str().c_str();
      if (name == "jd_et") {
         sscanf(value,"%lf",&(input.jd_et));
      }
      else if (name == "lon") {
         sscanf(value,"%lf",&(input.lon));
      }
      else if (name == "lat") {
         sscanf(value,"%lf",&(input.lat));
      }
      else if (name == "from") {
         sscanf(value,"%lf",&(input.from));
      }
      else if (name == "to") {
         sscanf(value,"%lf",&(input.to));
      }
      else if (name == "filter") {
         sscanf(value,"%d",&(input.filter));
      }
    }

}

int main( int argc, char** argv) {

  // Some defaults
  t_input input = { 2451545.,7.,52.,2452500.,2452501.,0b11111 };

  if (argc>=2) {
    // For tests
    parse_query_string(argv[1],input);
  } else {
    // In CGI mode
    const char* query_string = getenv("QUERY_STRING");
    if (query_string) {
      parse_query_string(query_string,input);
    }
  }

  cout << "HTTP/1.1 200 OK\nContent-Type:application/json;charset=utf-8\n\n";

  Horoscope h(input.jd_et,input.lon,input.lat);

// Default filter: 0b11111 = Down to and including ⚹
  Aspectarium a = h.getAspectarium( input.filter );
  multimap<double,Transit> t = transits(input.from,input.to,a);

// Output
  cout << "{ \"data\":" ;
  print(cout,t);
  cout << "}" << endl;

// Teardown - becomes relevant when code is shifted into shared library
  teardown( );

  return 0;

}