Interplay of Range and Pattern

2009-03-11
2013-04-09
  • The interplay of Range and Pattern is illustrated by the following prototypical example. Again, the file is essentially self-contained and only requires the Swiss Ephemeris library.

    In this simple example, an interval range is declared to step from 2000 till 2010 in 1-day intervals. The time for the check should be 6h ET, so we have a chance for the Sun to be in house I. We check the famous "Troinski" rule "Sun in house I, Mars in house XII, Saturn in X":

    // Example pattern: Troinski
    bool troinski( t_horo *h) {
     
      return ( ( house_pos( SE_SATURN, h) == 10 ) &&
               ( house_pos( SE_MARS, h)   == 12 ) &&
               ( house_pos( SE_SUN, h)    == 1  ) );
            
      }    

    The rule is attached as a C function pointer, returning true or false, to a Pattern object. We choose C function pointers here, because one option would be to embed Fabrice Bellards Tiny C Compiler TCC for an on-the-fly definition and execution of C functions at runtime.

    Observe the loop structure in main():

        while (r->hasMore()) {           // Iterate...
           r->getNextHoroscope( &h );    // Get next horoscope data (jd_et, lon, lat)
           if (p->applies(&h)) {         // Check the pattern
             swe_revjul( h.jd_et, 1, &year, &month, &day, &et );  // Calendar date                         
             printf( "%d.%d.%d (%lf)\n", day, month, year, h.jd_et );  // Output
             }    
           }  

    The output shows two intervals for which the condition is satisfied. 124 millisec execution is for a non-optimized code:

    -----------------------------------------------------
    18.7.2000 (2451743.750000)
    ...
    19.8.2000 (2451775.750000)
    19.8.2002 (2452505.750000)
    (...)
    4.9.2002 (2452521.750000)
    Execution time:    0.124 sec
    Press any key to continue . . .
    -----------------------------------------------------

    Source code of proto.cpp:

    #define USE_DLL
    #include "\sweph\src\swephexp.h"
    #include <time.h>
    #include <string>

    using namespace std;

    // -------------------------------------------------------------------
    // Example showing the interplay of Range and Proto
    // -------------------------------------------------------------------

    // Raw data of a horoscope = the triple (jd_et, lon, lat)
    // will eventually be part of a class Horoscope
    typedef struct {
      double jd_et;
      double lon;
      double lat;       
      } t_horo;

    // -------------------------------------------------------------------
    // Pattern
    // Base class checks for a score
    // -------------------------------------------------------------------
    class Pattern {
      public:        
        Pattern(bool (*)( t_horo *h ));                  
        bool applies(t_horo *h);
      protected:   
        bool (*rule)( t_horo *h );        
    };     

    // -------------------------------------------------------------------
    // Range
    // An arbitrary range of times or horoscopes
    // -------------------------------------------------------------------
    class Range {
      public:
        virtual double getNext() = 0;
        virtual bool getNextHoroscope(t_horo* h);
        virtual bool hasMore() = 0;
        virtual void set_place(double lon, double lat);
        double get_jd_start();
        double get_jd_end();
      protected:
        double jd_start, jd_end;                        
        bool place_provided;
        double lon, lat;  // Longitude and latitude if place is provided
      };     

    // -------------------------------------------------------------------
    // IntervalRange
    // Stepping through an interval with fixed step size
    // -------------------------------------------------------------------
    class IntervalRange : public Range {
      public:
        IntervalRange(double jd_start, double jd_end, double incr);              
        virtual double getNext();
        virtual bool hasMore();
      protected:
        double incr;                 
        double jd_cur;
      };     

    //-------------------------------------------------------------------- 
    //  --- Implementations 
    //-------------------------------------------------------------------- 

    Pattern::Pattern( bool (*rule) (t_horo *h) ){
      this->rule = rule;                 
      }
                     
    bool Pattern::applies( t_horo *h ) {
      if (this->rule == NULL)
        return false;
      else
        return (*rule)(h);    
      }    

    double Range::get_jd_start() {
      return this->jd_start;
      }    

    double Range::get_jd_end() {
      return this->jd_end;
      }    
     
    void Range::set_place(double lon, double lat) {
      this->lon = lon;
      this->lat = lat;
      this->place_provided = true;
      }      

    IntervalRange::IntervalRange(double jd_start, double jd_end, double incr) {
      this->jd_start = jd_cur = jd_start;                                   
      this->jd_end   = jd_end;                                   
      this->incr     = incr;                                   
      }                                                                          

    double IntervalRange::getNext() {
      double jd = jd_cur;
      jd_cur += incr;
      return jd;    
      }      

    bool IntervalRange::hasMore() {
      return (jd_cur <= jd_end);
      }      
     
    bool Range::getNextHoroscope(t_horo *h) {
      h->jd_et = getNext();
      if (place_provided) {
        h->lon = this->lon;
        h->lat = this->lat;
        return true;   
        }          
      else return false;               
      }     
     
    // House position of a planet 
    int house_pos( int ipl, t_horo *h) {
      double xx[6], armc, eps, jd_ut;
      char serr[255];
      int iflag = 0, iret;
      jd_ut = h->jd_et - swe_deltat(h->jd_et);
      armc = swe_sidtime(jd_ut) * 15;  
      swe_calc(h->jd_et, SE_ECL_NUT, 0, xx, serr);
      eps = xx[0];
      iret = swe_calc( h->jd_et, ipl, iflag, xx, serr );
      xx[1] = 0.;  // ignore ecliptical latitude
      return (int) swe_house_pos( armc, h->lat, eps, 'P', xx, serr);
      }       
      

    // Example pattern: Troinski
    bool troinski( t_horo *h) {
     
      return ( ( house_pos( SE_SATURN, h) == 10 ) &&
               ( house_pos( SE_MARS, h)   == 12 ) &&
               ( house_pos( SE_SUN, h)    == 1  ) );
            
      }    

    int main(int argc, char *argv[])
    {
       
        int year1 = 2000,
            year2 = 2010,
            year, month, day;
        double et;       
        t_horo h = { 0., 0., 50. };     
        Range *r;
        Pattern *p = new Pattern( & troinski );
       
        double starttime = (double)clock(), endtime;

        r = new IntervalRange(
          swe_julday(year1,1,1,6.,1),
          swe_julday(year2,1,1,0.,1),
          1);   
         
         
        while (r->hasMore()) {           // Iterate...
           r->getNextHoroscope( &h );    // Get next horoscope data (jd_et, lon, lat)
           if (p->applies(&h)) {         // Check the pattern
             swe_revjul( h.jd_et, 1, &year, &month, &day, &et );  // Calendar date                         
             printf( "%d.%d.%d (%lf)\n", day, month, year, h.jd_et );  // Output
             }    
           }  

        endtime = (double) clock();
        printf("Execution time: %8.3lf sec\n",(endtime - starttime)/CLOCKS_PER_SEC);

        system("PAUSE");
        return EXIT_SUCCESS;
    }

     
    • To complete the example, I added the ingress determination functionality (for arbitrary planets, although only used for the Sun in this example).

      #define USE_DLL
      #include "\sweph\src\swephexp.h"
      #include <time.h>
      #include <string>

      // -------------------------------------------------------------------
      // ad hoc functions for default determination of the gregflag
      // -------------------------------------------------------------------
      #define gregflag_from_year(x) (x > 1582 ? 1 : 0)
      #define gregflag_from_jd(x) (x > 2298883.5 ? 1 : 0)
      #define ut_from_jd_et(x)((x-swe_deltat(x)+0.5- floor(x-swe_deltat(x)+0.5))*24 )
      //#define TRACE

      using namespace std;

      // -------------------------------------------------------------------
      // Example showing the interplay of Range and Proto
      // -------------------------------------------------------------------

      // Raw data of a horoscope = the triple (jd_et, lon, lat)
      // will eventually be part of a class Horoscope
      typedef struct {
        double jd_et;
        double lon;
        double lat;
        } t_horo;

      // -------------------------------------------------------------------
      // Pattern
      // Base class checks for a score
      // -------------------------------------------------------------------
      class Pattern {
        public:        
          Pattern(bool (*)( t_horo *h ));                  
          bool applies(t_horo *h);
        protected:   
          bool (*rule)( t_horo *h );        
      };     

      // -------------------------------------------------------------------
      // Range
      // An arbitrary range of times or horoscopes
      // -------------------------------------------------------------------
      class Range {
        public:
          virtual double getNext() = 0;
          virtual bool getNextHoroscope(t_horo* h);
          virtual bool hasMore() = 0;
          virtual void set_place(double lon, double lat);
          double get_jd_start();
          double get_jd_end();
        protected:
          double jd_start, jd_end;                        
          bool place_provided;
          double lon, lat;  // Longitude and latitude if place is provided
        };     

      // -------------------------------------------------------------------
      // Valuation (abstract class)
      // Determine zeroes of a function f implemented in a subclass
      // -------------------------------------------------------------------
      class Valuation {
        public:   
      // Target function, has to be (re)defined:
          virtual double f(double x) = 0;
      // If some intervals can be excluded with a cheap check,
      // the scan will become more efficient:  
          virtual bool skip(double x0, double x1, double y0, double y1) { return false; };
      // The algorithm itself (Regula falsi):   
          bool getZero( double x0, double x1, double *x  );   
        private:
      // precision in time for events the millionth part of a day
      // sufficient for precision to a second
          static const double PREC = 1.E-6;
      // precision in longitude (for the result)
          static const double LPREC = 1.E-10;
      // maximal number of iterations (to prevent endless loops)
          static const int MAXITER = 100;             
        };

      // -------------------------------------------------------------------
      // IntervalRange
      // Stepping through an interval with fixed step size
      // -------------------------------------------------------------------
      class IntervalRange : public Range {
        public:
          IntervalRange(double jd_start, double jd_end, double incr);              
          virtual double getNext();
          virtual bool hasMore();
        protected:
          double incr;                 
          double jd_cur;
        };     
       

      // -------------------------------------------------------------------
      // Ingresses (planet entering new zodiac sign)
      // -------------------------------------------------------------------
      class IngressRange : public Range, Valuation {
        public:
          IngressRange( int ipl, double jd_start, double jd_end);            
          virtual double getNext();
          virtual bool hasMore();
          virtual double f(double x);
        protected:
          int ipl;           
          double jd_cur;            
          double incr;
          double reductionValue; // used in valuation function f
          bool inRange;
        };     

      // -------------------------------------------------------------------
      // Cardinal Ingresses (planet entering new zodiac sign)
      // -------------------------------------------------------------------
      class CardinalIngressRange : public IngressRange {
        public:   
          CardinalIngressRange(int ipl, double jd_start, double jd_end);
        };     

      // -------------------------------------------------------------------
      // Aries Ingresses (planet entering Aries)
      // -------------------------------------------------------------------
      class AriesIngressRange : public IngressRange {
        public:   
          AriesIngressRange(int ipl, double jd_start, double jd_end);
          virtual double f(double x);
          virtual bool skip(double x0, double x1, double y0, double y1);
        };     

      //-------------------------------------------------------------------- 
      //  --- Implementations 
      //-------------------------------------------------------------------- 

      Pattern::Pattern( bool (*rule) (t_horo *h) ){
        this->rule = rule;                 
        }
                       
      bool Pattern::applies( t_horo *h ) {
        if (this->rule == NULL)
          return false;
        else
          return (*rule)(h);    
        }    

      double Range::get_jd_start() {
        return this->jd_start;
        }    

      double Range::get_jd_end() {
        return this->jd_end;
        }    
       
      void Range::set_place(double lon, double lat) {
        this->lon = lon;
        this->lat = lat;
        this->place_provided = true;
        }      

      IntervalRange::IntervalRange(double jd_start, double jd_end, double incr) {
        this->jd_start = jd_cur = jd_start;                                   
        this->jd_end   = jd_end;                                   
        this->incr     = incr;                                   
        }                                                                          

      double IntervalRange::getNext() {
        double jd = jd_cur;
        jd_cur += incr;
        return jd;    
        }      

      bool IntervalRange::hasMore() {
        return (jd_cur <= jd_end);
        }      
       
      bool Range::getNextHoroscope(t_horo *h) {
        h->jd_et = getNext();
        if (place_provided) {
          h->lon = this->lon;
          h->lat = this->lat;
          return true;   
          }          
        else return false;               
        }     
       
      //-------------------------------------------------------------------- 
      // Ingresses
      //-------------------------------------------------------------------- 
      IngressRange::IngressRange(int ipl, double jd_start, double jd_end)
      {
        this->ipl      = ipl;                              
        this->jd_start = jd_cur = jd_start;                                   
        this->jd_end   = jd_end;
      // incr must be smaller than 15/MaxSpeed(°)(ipl) days 
        if (ipl == SE_SUN  || ipl >= SE_MARS) {
          this->incr = 15.;         
          }        
        else if (ipl == SE_MOON) {
          this->incr = 0.9;      
          }
        else {
          this->incr = 3.;   
          }                     
        this->reductionValue = 15.;  // used in valuation function
        this->place_provided = false;                              
        }
                                    
      // Valuation for an ingress
      double IngressRange::f( double jd ) {
        double xx[6],x,y;
        int iflag = 0, iret;
        char serr[255];
        iret = swe_calc( jd, ipl, iflag, xx, serr);    
        x = xx[0] / reductionValue;
        y =  1. - fabs( x - 4*floor( (x+1.)/4 ) - 1.);
        return y;
        }      
       
      bool IngressRange::hasMore() {
        double x1,x;
        for (int i=0;; i++) {
          jd_cur += incr;
          if ( jd_cur > jd_end ) return (inRange = false);
          x1 = jd_cur + incr;
          if ( x1 > jd_end ) x1 = jd_end;
          inRange = getZero( jd_cur,x1, &x);
          if (inRange) {
            jd_cur = x;           
            return true;
            }
          }
        return (inRange = false);              
        }    
       
      double IngressRange::getNext() {
        return inRange ? jd_cur : 0;    
        }        
       
      CardinalIngressRange::CardinalIngressRange(int ipl, double jd_start, double jd_end) :
            IngressRange(ipl,jd_start,jd_end) {
        this->incr     = 3 * incr; 
        this->reductionValue = 45.;  // used in valuation function
        }                                         

      AriesIngressRange::AriesIngressRange(int ipl, double jd_start, double jd_end) :
            IngressRange(ipl,jd_start,jd_end) {
        this->incr           = 50.; 
        this->reductionValue = 180.;  // used in valuation function
        }   
       
      double AriesIngressRange::f(double jd) {
        double xx[6],y;
        int iflag = 0, iret;
        char serr[255];
        iret = swe_calc( jd, ipl, iflag, xx, serr);    
        y = xx[0] / 360.;
        return y - floor( y + 0.5 );    
        }      
       
      bool AriesIngressRange::skip(double x0, double x1, double y0, double y1) {
        return (fabs(y0) > .1) && (fabs(y1) > .1);       
        }

      //-------------------------------------------------------------------- 
      // Determine a zero of a function
      //-------------------------------------------------------------------- 
      bool Valuation::getZero( double _x0, double _x1, double *x ) {
                            
        double x0 = _x0,
               x1 = _x1,
               y0,y1,y;
         int iter = 0;    
        
      // Check the endpoints first  
         if (fabs(y0 = f(x0)) < LPREC) {
           *x = x0;
           return true;
           }
         if (fabs(y1 = f(x1)) < LPREC)  {
           *x = x1;
           return true;
           }

      // Cheap check...
         if (skip(x0,x1,y0,y1)) return false;
          
      // Provided the interval is sufficiently small (this falls into
      // responsibility of the caller), a zero requires a sign change
      // of the function values at the endpoints
          if ( ( (y0 > 0) && (y1 > 0) ) ||
               ( (y0 < 0) && (y1 < 0) ) )
            
             return false;
              
      // Use Regula Falsi to get the zero
          do {
            if (++iter > MAXITER) {
              fprintf(stderr,"Iterationsfehler\n");
              return false;
              }

            *x = x0 - y0*(x1-x0)/(y1-y0);
            y = f( *x );

      #ifdef TRACE
              printf(" d0:%12.6f; d1:%12.6f; y:%12.6f\n",y0,y1,y);
              printf(" x0:%12.6f; x1:%12.6f; x:%12.6f\n",x0,x1,*x);
      #endif

            if (fabs(y) < LPREC) break;
            if ((y<0)&&(y0>0) || (y>0)&&(y0<0)) {
              y1  = y;
              x1  = *x;
              }
            else {
              y0 = y;
              x0 = *x; 
              } 
            }
          while (fabs(x1-x0)>PREC);
          return true;
          }     

                                                       
                                                        
      // House position of a planet 
      int house_pos( int ipl, t_horo *h) {
        double xx[6], armc, eps, jd_ut;
        char serr[255];
        int iflag = 0, iret;
        jd_ut = h->jd_et - swe_deltat(h->jd_et);
        armc = swe_sidtime(jd_ut) * 15 + h->lon;  
        swe_calc(h->jd_et, SE_ECL_NUT, 0, xx, serr);
        eps = xx[0];
        iret = swe_calc(  h->jd_et, ipl, iflag, xx, serr );
        xx[1] = 0.;  // ignore ecliptical
        return (int) swe_house_pos( armc, h->lat, eps, 'P', xx, serr);
        }       
        

      // Example pattern: Troinski
      bool troinski( t_horo *h) {
         
        return ( ( house_pos( SE_SATURN, h) == 10 ) &&
                 ( house_pos( SE_MARS, h)   == 12 ) &&
                 ( house_pos( SE_SUN, h)    == 1  ) );
              
        }    

      // --- Build a range from the data given in the command line
      void parse_command_line( int argc, char* argv[], Range **r) {

        char ingress_type;
        int year1, year2, ipl;
        double jd_start, jd_end;
       
          if (argc < 4) {
            printf("Usage: proto [aci] <planet> <year_from> <year_to>\n");
            exit(-1);
            }
          else {
            sscanf( argv[1], "%c", &ingress_type );
            sscanf( argv[2], "%d", &ipl );
            sscanf( argv[3], "%d", &year1 );
            jd_start = swe_julday(year1,1,1,0.,gregflag_from_year(year1));
            sscanf( argv[4], "%d", &year2 );
            jd_end =  swe_julday(year2,1,1,0.,gregflag_from_year(year1));     
            }  
         
          switch( ingress_type ) {
            case 'i' :            
              *r = new IngressRange( ipl, jd_start, jd_end );
              break;
            case 'c' :            
              *r = new CardinalIngressRange( ipl, jd_start, jd_end );
              break;
            case 'a' :            
              *r = new AriesIngressRange( ipl, jd_start, jd_end );
              break;
            }       

      }    

      //-------------------------------------------------------------------- 
      int main(int argc, char *argv[])
      {
          int year, month, day;
          double et, starttime = (double)clock(), endtime;
          t_horo h = { 0., 12.5, 41.9 };     
          Range *r;
         
      // Search for the Troinski pattern "Sun in I, Mars in XII, Saturn in X"   
          Pattern *p = new Pattern( & troinski );

      // Parse command line options
          parse_command_line( argc, argv, &r);
           
          while (r->hasMore()) {           // Iterate...
             r->getNextHoroscope( &h );    // Get next horoscope data (jd_et, lon, lat)
             if (p->applies(&h)) {         // Check the pattern
               swe_revjul( h.jd_et, gregflag_from_jd(h.jd_et), &year, &month, &day, &et );  // Calendar date                         
               printf( "%2d.%2d.%5d %6.2lfh UT\n", day, month, year, ut_from_jd_et(h.jd_et) );  // Output
               }    
             }  

          endtime = (double) clock();
          printf("Execution time: %8.3lf sec\n",(endtime - starttime)/CLOCKS_PER_SEC);

          system("PAUSE");
          return EXIT_SUCCESS;
      }