/* SWISSEPH 

  Heliacal risings and related calculations
  
  Author: Victor Reijs
  This program code is a translation of part of:
  Victor Reijs' software ARCHAEOCOSMO (archaeoastronomy and
  geodesy functions), 
  http://www.iol.ie/~geniet/eng/archaeocosmoprocedures.htm

  Translation from VB into C by Dieter Koch

  Problem reports can be sent to victor.reijs@gmail.com or dieter@astro.ch
  
  Copyright (c) Victor Reijs, 2008
  Copyright (C) 1997 - 2021 Astrodienst AG, Switzerland.  All rights reserved.

  License conditions
  ------------------

  This file is part of Swiss Ephemeris.

  Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND.  No author
  or distributor accepts any responsibility for the consequences of using it,
  or for whether it serves any particular purpose or works at all, unless he
  or she says so in writing.  

  Swiss Ephemeris is made available by its authors under a dual licensing
  system. The software developer, who uses any part of Swiss Ephemeris
  in his or her software, must choose between one of the two license models,
  which are
  a) GNU Affero General Public License (AGPL)
  b) Swiss Ephemeris Professional License

  The choice must be made before the software developer distributes software
  containing parts of Swiss Ephemeris to others, and before any public
  service using the developed software is activated.

  If the developer choses the AGPL software license, he or she must fulfill
  the conditions of that license, which includes the obligation to place his
  or her whole software project under the AGPL or a compatible license.
  See https://www.gnu.org/licenses/agpl-3.0.html

  If the developer choses the Swiss Ephemeris Professional license,
  he must follow the instructions as found in http://www.astro.com/swisseph/ 
  and purchase the Swiss Ephemeris Professional Edition from Astrodienst
  and sign the corresponding license contract.

  The License grants you the right to use, copy, modify and redistribute
  Swiss Ephemeris, but only under certain conditions described in the License.
  Among other things, the License requires that the copyright notices and
  this notice be preserved on all copies.

  Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl

  The authors of Swiss Ephemeris have no control or influence over any of
  the derived works, i.e. over software or services created by other
  programmers which use Swiss Ephemeris functions.

  The names of the authors or of the copyright holder (Astrodienst) must not
  be used for promoting any software, product or service which uses or contains
  the Swiss Ephemeris. This copyright notice is the ONLY place where the
  names of the authors can legally appear, except in cases where they have
  given special permission in writing.

  The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
  for promoting such software, products or services.
*/


#include "swephexp.h"
#include "sweph.h"
#include "swephlib.h"
#include <sys/stat.h>

#define PLSV   0 /*if Planet, Lunar and Stellar Visibility formula is needed PLSV=1*/
#define criticalangle   0.0 /*[deg]*/
#define BNIGHT   	1479.0 /*[nL]*/
#define BNIGHT_FACTOR   1.0
#define PI		M_PI
#define Min2Deg   (1.0 / 60.0)
#define SWEHEL_DEBUG  0
#define DONE  1
#define MaxTryHours   4 
#define TimeStepDefault	1
#define LocalMinStep	8

/* time constants */
#define Y2D   365.25 /*[Day]*/
#define D2Y   (1 / Y2D) /*[Year]*/
#define D2H   24.0 /*[Hour]*/
#define H2S   3600.0 /*[sec]*/
#define D2S   (D2H * H2S) /*[sec]*/
#define S2H   (1.0 / H2S) /*[Hour]*/
#define JC2D   36525.0 /*[Day]*/
#define M2S   60.0 /*[sec]*/

/* Determines which algorimths are used*/
#define REFR_SINCLAIR    0
#define REFR_BENNETTH    1
#define FormAstroRefrac   REFR_SINCLAIR /*for Astronomical refraction can be "bennetth" or "sinclair"*/
#define GravitySource   2 /*0=RGO, 1=Wikipedia,2=Exp. Suppl. 1992,3=van der Werf*/
#define REarthSource   1 /*0=RGO (constant), 1=WGS84 method*/

#define StartYear   1820 /*[year]*/
#define Average   1.80546834626888 /*[msec/cy]*/
#define Periodicy   1443.67123144531 /*[year]*/
#define Amplitude   3.75606495492684 /*[msec]*/
#define phase   0 /*[deg]*/
#define MAX_COUNT_SYNPER           5  /* search within 10 synodic periods */
#define MAX_COUNT_SYNPER_MAX 1000000  /* high, so there is not max count */
#define AvgRadiusMoon  (15.541 / 60) /* '[Deg] at 2007 CE or BCE*/

/* WGS84 ellipsoid constants
 * http://w3sli.wcape.gov.za/Surveys/Mapping/wgs84.htm*/
#define Ra   6378136.6   /*'[m]*/
#define Rb   6356752.314 /*'[m]*/

/* choices in Schaefer's model */
#define nL2erg 	(1.02E-15)
#define erg2nL 	(1 / nL2erg) /*erg2nL to nLambert*/
#define MoonDistance 	384410.4978 /*[km]*/
#define scaleHwater 	3000.0 /*[m] Ricchiazzi [1997] 8200 Schaefer [2000]*/
#define scaleHrayleigh 	8515.0 /*[m] Su [2003] 8200 Schaefer [2000]*/
#define scaleHaerosol 	3745.0 /*m Su [2003] 1500 Schaefer [2000]*/
#define scaleHozone 	20000.0 /*[m] Schaefer [2000]*/
#define astr2tau 	0.921034037197618  /*LN(10 ^ 0.4)*/
#define tau2astr 	1 / astr2tau

/* meteorological constants*/
#define C2K   273.15 /*[K]*/
#define DELTA   18.36
#define TempNulDiff   0.000001
#define PressRef   1000 /*[mbar]*/
#define MD   28.964 /*[kg] Mol weight of dry air van der Werf*/
#define MW   18.016 /*[kg] Mol weight of water vapor*/
#define GCR   8314.472 /*[L/kmol/K] van der Werf*/
#define LapseSA   0.0065 /*[K/m] standard atmosphere*/
#define LapseDA   0.0098 /*[K/m] dry adiabatic*/

/* lowest apparent altitude to provide*/
#define LowestAppAlt   -3.5 /*[Deg]*/

/*optimization delta*/
#define epsilon   0.001
/* for Airmass usage*/
#define staticAirmass   0 /* use staticAirmass=1 instead depending on difference k's*/

/* optic stuff */
#define GOpticMag   1 /*telescope magnification*/
#define GOpticTrans   0.8 /*telescope transmission*/
#define GBinocular   1 /*1-binocular 0=monocular*/
#define GOpticDia   50 /*telescope diameter [mm]*/

static double mymin(double a, double b)
{
  if (a <= b)
    return a;
  return b;
}

static double mymax(double a, double b)
{
  if (a >= b)
    return a;
  return b;
}

/*###################################################################*/
static double Tanh(double x)
{
  return (exp(x) - exp(-x)) / (exp(x) + exp(-x));
}

/*
' B [nL]
' SN [-]
' CVA [deg]
*/
static double CVA(double B, double SN, int32 helflag)
{
  /*Schaefer, Astronomy and the limits of vision, Archaeoastronomy, 1993*/
  AS_BOOL is_scotopic = FALSE;
  //if (B < BNIGHT)
  if (B < 1394)  /* use this value for BNIGHT to make the function continous */
    is_scotopic = TRUE;
  if (helflag & SE_HELFLAG_VISLIM_PHOTOPIC) 
    is_scotopic = FALSE;
  if (helflag & SE_HELFLAG_VISLIM_SCOTOPIC) 
    is_scotopic = TRUE;
  if (is_scotopic)
    return mymin(900, 380 / SN * pow(10, (0.3 * pow(B, (-0.29))))) / 60.0 / 60.0;
  else
    return (40.0 / SN) * pow(10, (8.28 * pow(B, (-0.29)))) / 60.0 / 60.0;
}

/*
' age [year]
' B [nL]
' PupilDia [mm]
*/
static double PupilDia(double Age, double B)
{
  /* age dependancy from Garstang [2000]*/
  return (0.534 - 0.00211 * Age - (0.236 - 0.00127 * Age) * Tanh(0.4 * log(B) / log(10) - 2.2)) * 10;
}

/*
'Input
' Bback [nL]
' kX [-]
' Binocular [-]
' OpticMag [-]
' OpticDia [mm]
' OpticTrans [-]
' JDNDaysUT [JDN]
' Age [Year]
' SN [-]
' ObjectName
' TypeFactor [0=itensity factor 1=background factor]
'Output
' OpticFactor [-]
*/
static double OpticFactor(double Bback, double kX, double *dobs, double JDNDaysUT, char *ObjectName, int TypeFactor, int helflag)
{
  double Pst, CIb, CIi, ObjectSize, Fb, Fe, Fsc, Fci, Fcb, Ft, Fp, Fa, Fr, Fm;
  double Age = dobs[0];
  double SN = dobs[1], SNi;
  double Binocular = dobs[2];
  double OpticMag = dobs[3];
  double OpticDia = dobs[4];
  double OpticTrans = dobs[5];
  AS_BOOL is_scotopic = FALSE;
  JDNDaysUT += 0.0; /* currently not used, statement prevents compiler warning */
  SNi = SN;
  if (SNi <= 0.00000001) SNi = 0.00000001;
  /* 23 jaar as standard from Garstang*/
  Pst = PupilDia(23, Bback);
  if (OpticMag == 1) { /*OpticMagn=1 means using eye*/
     OpticTrans = 1;
     OpticDia = Pst;
  }
#if 0 /*is done in default_heliacal_parameters()*/
  if (OpticMag == 0) { /*OpticMagn=0 (undefined) using eye*/
    OpticTrans = 1;
    OpticDia = Pst;
    Binocular = 1;
    OpticMag = 1;
  }
#endif
  /* Schaefer, Astronomy and the limits of vision, Archaeoastronomy, 1993*/
  CIb = 0.7; /* color of background (from Ben Sugerman)*/
  CIi = 0.5; /* Color index for white (from Ben Sugerman), should be function of ObjectName*/
  ObjectSize = 0;
  if (strcmp(ObjectName, "moon") == 0) {
    /*ObjectSize and CI needs to be determined (depending on JDNDaysUT)*/
    ;
  }
  Fb = 1;
  if (Binocular == 0) Fb = 1.41;
  //if (Bback < BNIGHT)
  if (Bback < 1645) /* use this value for BNIGHT to make the function continuous */
    is_scotopic = TRUE;
  if (helflag & SE_HELFLAG_VISLIM_PHOTOPIC) 
    is_scotopic = FALSE;
  if (helflag & SE_HELFLAG_VISLIM_SCOTOPIC) 
    is_scotopic = TRUE;
  if (is_scotopic) {
    Fe = pow(10, (0.48 * kX));
    Fsc = mymin(1, (1 - pow(Pst / 124.4, 4)) / (1 - pow((OpticDia / OpticMag / 124.4), 4)));
    Fci = pow(10, (-0.4 * (1 - CIi / 2.0)));
    Fcb = pow(10, (-0.4 * (1 - CIb / 2.0)));
  } else {
     Fe = pow(10, (0.4 * kX));
     Fsc = mymin(1, pow((OpticDia / OpticMag / Pst), 2) * (1 - exp(-pow((Pst / 6.2), 2))) / (1 - exp(-pow((OpticDia / OpticMag / 6.2), 2))));
     Fci = 1;
     Fcb = 1;
  }
  Ft = 1 / OpticTrans;
  Fp = mymax(1, pow((Pst / (OpticMag * PupilDia(Age, Bback))), 2));
  Fa = pow((Pst / OpticDia), 2);
  Fr = (1 + 0.03 * pow((OpticMag * ObjectSize / CVA(Bback, SNi, helflag)), 2)) / pow(SNi, 2);
  Fm = pow(OpticMag, 2);
#if SWEHEL_DEBUG
  fprintf(stderr, "Pst=%f\n", Pst);
  fprintf(stderr, "Fb =%f\n", Fb);
  fprintf(stderr, "Fe =%f\n", Fe);
  fprintf(stderr, "Ft =%f\n", Ft);
  fprintf(stderr, "Fp =%f\n", Fp);
  fprintf(stderr, "Fa =%f\n", Fa);
  fprintf(stderr, "Fm =%f\n", Fm);
  fprintf(stderr, "Fsc=%f\n", Fsc);
  fprintf(stderr, "Fci=%f\n", Fci);
  fprintf(stderr, "Fcb=%f\n", Fcb);
  fprintf(stderr, "Fr =%f\n", Fr );
#endif
  if (TypeFactor == 0) 
    return Fb * Fe * Ft * Fp * Fa * Fr * Fsc * Fci;
  else
    return Fb * Ft * Fp * Fa * Fm * Fsc * Fcb;
}

/*###################################################################
*/
static int32 DeterObject(char *ObjectName)
{
  char s[AS_MAXCH];
  char *sp;
  int32 ipl;
  strcpy(s, ObjectName);
  for (sp = s; *sp != '\0'; sp++)
    *sp = tolower(*sp);
  if (strncmp(s, "sun", 3) == 0)
    return SE_SUN;
  if (strncmp(s, "venus", 5) == 0)
    return SE_VENUS;
  if (strncmp(s, "mars", 4) == 0)
    return SE_MARS;
  if (strncmp(s, "mercur", 6) == 0)
    return SE_MERCURY;
  if (strncmp(s, "jupiter", 7) == 0)
    return SE_JUPITER;
  if (strncmp(s, "saturn", 6) == 0)
    return SE_SATURN;
  if (strncmp(s, "uranus", 6) == 0)
    return SE_URANUS;
  if (strncmp(s, "neptun", 6) == 0)
    return SE_NEPTUNE;
  if (strncmp(s, "moon", 4) == 0)
    return SE_MOON;
  if ((ipl = atoi(s)) > 0) {
    ipl += SE_AST_OFFSET;
    return ipl;
  }
  return -1;
}

#if 0
int32 call_swe_calc(double tjd, int32 ipl, int32 iflag, double *x, char *serr) 
{
  int32 retval = OK, ipli, i;
  double dtjd;
  static TLS double tjdsv[3];
  static TLS double xsv[3][6];
  static TLS int32 iflagsv[3];
  ipli = ipl;
  if (ipli > SE_MOON) 
    ipli = 2;
  dtjd = tjd - tjdsv[ipli];
  if (tjdsv[ipli] != 0 && iflag == iflagsv[ipli] && fabs(dtjd) < 5.0 / 1440.0) {
    for (i = 0; i < 3; i++) 
      x[i] = xsv[ipli][i] + dtjd * xsv[ipli][i+3];
    for (i = 3; i < 6; i++) 
      x[i] = xsv[ipli][i];
  } else {
    retval = swe_calc(tjd, ipl, iflag, x, serr);
    tjdsv[ipli] = tjd;
    iflagsv[ipli] = iflag;
    for (i = 0; i < 6; i++) 
      xsv[ipli][i] = x[i];
  }
  return retval;
}
#endif

/* avoids problems with star name string that may be overwritten by 
   swe_fixstar() */
static int32 call_swe_fixstar(char *star, double tjd, int32 iflag, double *xx, char *serr)
{
  int32 retval;
  char star2[AS_MAXCH];
  strcpy(star2, star);
  retval =  swe_fixstar(star2, tjd, iflag, xx, serr);
  return retval;
}

/* avoids problems with star name string that may be overwritten by 
   swe_fixstar_mag() */
static int32 call_swe_fixstar_mag(char *star, double *mag, char *serr)
{
  int32 retval;
  char star2[AS_MAXCH];
  static TLS double dmag;
  static TLS char star_save[AS_MAXCH];
  if (strcmp(star, star_save) == 0) {
    *mag = dmag;
    return OK;
  }
  strcpy(star_save, star);
  strcpy(star2, star);
  retval = swe_fixstar_mag(star2, &dmag, serr);
  *mag = dmag;
  return retval;
}

/* avoids problems with star name string that may be overwritten by 
   swe_fixstar() */
static int32 call_swe_rise_trans(double tjd, int32 ipl, char *star, int32 helflag, int32 eventtype, double *dgeo, double atpress, double attemp, double *tret, char *serr)
{
  int32 retval;
  int32 iflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  char star2[AS_MAXCH];
  strcpy(star2, star);
  retval = swe_rise_trans(tjd, ipl, star2, iflag, eventtype, dgeo, atpress, attemp, tret, serr);
  return retval;
}

/* 
 * Written by Dieter Koch:
 * Fast function for risings and settings of planets, can be used instead of 
 * swe_rise_trans(), which is much slower.
 * For circumpolar and near-circumpolar planets use swe_rise_trans(), or 
 * generally use it for geographical latitudes higher than 58N/S.
 * For fixed stars, swe_rise_trans() is fast enough.
 */
static int32 calc_rise_and_set(double tjd_start, int32 ipl, double *dgeo, double *datm, int32 eventflag, int32 helflag, double *trise, char *serr)
{
  int retc = OK, i;
  double sda, xs[6], xx[6], xaz[6], xaz2[6], dfac = 1/365.25;
  double rdi, rh;
  double tjd0 = tjd_start, tjdrise;
  double tjdnoon = (int) tjd0 - dgeo[0] / 15.0 / 24.0;
  int32 iflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  int32 epheflag = iflag;
  iflag |= SEFLG_EQUATORIAL;
  if (!(helflag & SE_HELFLAG_HIGH_PRECISION)) 
    iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
  if (swe_calc_ut(tjd0, SE_SUN, iflag, xs, serr) == 0) {
    if (serr != NULL)
      strcpy(serr, "error in calc_rise_and_set(): calc(sun) failed ");
    return ERR;
  }
  if (swe_calc_ut(tjd0, ipl, iflag, xx, serr) == 0) {
    if (serr != NULL)
      strcpy(serr, "error in calc_rise_and_set(): calc(sun) failed ");
    return ERR;
  }
  tjdnoon -= swe_degnorm(xs[0] - xx[0])/360.0 + 0;
  /* is planet above horizon or below? */
  swe_azalt(tjd0, SE_EQU2HOR, dgeo, datm[0], datm[1], xx, xaz);
if (eventflag & SE_CALC_RISE) {
  if (xaz[2] > 0) {
    while (tjdnoon - tjd0 < 0.5) {/*printf("e");*/tjdnoon += 1;}
    while (tjdnoon - tjd0 > 1.5) {/*printf("f");*/tjdnoon -= 1;}
  } else {
    while (tjdnoon - tjd0 < 0.0) {/*printf("g");*/tjdnoon += 1;}
    while (tjdnoon - tjd0 > 1.0) {/*printf("h");*/tjdnoon -= 1;}
  }
} else {
  if (xaz[2] > 0) {
    while (tjd0 - tjdnoon > 0.5) {/*printf("a");*/ tjdnoon += 1;}
    while (tjd0 - tjdnoon < -0.5) {/*printf("b");*/ tjdnoon -= 1;}
  } else {
    while (tjd0 - tjdnoon > 0.0) {/*printf("c");*/ tjdnoon += 1;}
    while (tjd0 - tjdnoon < -1.0) {/*printf("d");*/ tjdnoon -= 1;}
  }
}
  /* position of planet */
  if (swe_calc_ut(tjdnoon, ipl, iflag, xx, serr) == ERR) {
    if (serr != NULL)
      strcpy(serr, "error in calc_rise_and_set(): calc(sun) failed ");
    return ERR;
  }
  /* apparent radius of solar disk (ignoring refraction) */
  rdi = 0;
  if (ipl == SE_SUN) 
    rdi = asin(696000000.0 / 1.49597870691e+11 / xx[2]) / DEGTORAD;
  else if (ipl == SE_MOON)
    rdi = asin(1737000.0 / 1.49597870691e+11 / xx[2]) / DEGTORAD;
  if (eventflag & SE_BIT_DISC_CENTER)
    rdi = 0;
  /* true altitude of sun, when it appears at the horizon */
  /* refraction for a body visible at the horizon at 0m above sea,
   * atmospheric temperature 10 deg C, atmospheric pressure 1013.25 is 34.5 arcmin*/
  rh = -(34.5 / 60.0 + rdi);
  /* semidiurnal arc of sun */
  sda = acos(-tan(dgeo[1] * DEGTORAD) * tan(xx[1] * DEGTORAD)) * RADTODEG;
  /* rough rising and setting times */
  if (eventflag & SE_CALC_RISE)
    tjdrise = tjdnoon - sda / 360.0;
  else
    tjdrise = tjdnoon + sda / 360.0;
  /*ph->tset = tjd_start + sda / 360.0;*/
  /* now calculate more accurate rising and setting times.
   * use vertical speed in order to determine crossing of the horizon  
   * refraction of 34' and solar disk diameter of 16' = 50' = 0.84 deg */
  iflag = epheflag|SEFLG_SPEED|SEFLG_EQUATORIAL;
  if (ipl == SE_MOON)
    iflag |= SEFLG_TOPOCTR;
  if (!(helflag & SE_HELFLAG_HIGH_PRECISION)) 
    iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
  for (i = 0; i < 2; i++) {
    if (swe_calc_ut(tjdrise, ipl, iflag, xx, serr) == ERR) {
      /*fprintf(stderr, "hev4 tjd=%f, ipl=%d, iflag=%d\n", tjdrise, ipl, iflag);*/
      return ERR;
    }
    swe_azalt(tjdrise, SE_EQU2HOR, dgeo, datm[0], datm[1], xx, xaz);
    xx[0] -= xx[3] * dfac; 
    xx[1] -= xx[4] * dfac;
    swe_azalt(tjdrise - dfac, SE_EQU2HOR, dgeo, datm[0], datm[1], xx, xaz2);
    tjdrise -= (xaz[1] - rh) / (xaz[1] - xaz2[1]) *  dfac;
    /*fprintf(stderr, "%f\n", ph->trise);*/
  }
  *trise = tjdrise;
  return retc;
}

static int32 my_rise_trans(double tjd, int32 ipl, char* starname, int32 eventtype, int32 helflag, double *dgeo, double *datm, double *tret, char *serr)
{
  int retc = OK;
  if (starname != NULL && *starname != '\0')
    ipl = DeterObject(starname);
  /* for non-circumpolar planets we can use a faster algorithm */
  /*if (!(helflag & SE_HELFLAG_HIGH_PRECISION) && ipl != -1 && fabs(dgeo[1]) < 58) {*/
  if (ipl != -1 && fabs(dgeo[1]) < 63) {
    retc = calc_rise_and_set(tjd, ipl, dgeo, datm, eventtype, helflag, tret, serr);
  /* for stars and circumpolar planets we use a rigorous algorithm */
  } else {
    retc = call_swe_rise_trans(tjd, ipl, starname, helflag, eventtype, dgeo, datm[0], datm[1], tret, serr);
  }
/*  printf("%f, %f\n", tjd, *tret);*/
  return retc;
}

/*###################################################################
' JDNDaysUT [Days]
' dgeo [array: longitude, latitude, eye height above sea m]
' TempE [C]
' PresE [mbar]
' ObjectName (string)
' RSEvent (1=rise, 2=set,3=up transit,4=down transit)
' Rim [0=center,1=top]
' RiseSet [Day]
*/
static int32 RiseSet(double JDNDaysUT, double *dgeo, double *datm, char *ObjectName, int32 RSEvent, int32 helflag, int32 Rim, double *tret, char *serr)
{
  int32 eventtype = RSEvent, Planet, retval;
  if (Rim == 0)
    eventtype |= SE_BIT_DISC_CENTER;
  Planet = DeterObject(ObjectName);
  if (Planet != -1)
    retval = my_rise_trans(JDNDaysUT, Planet, "", eventtype, helflag, dgeo, datm, tret, serr);
  else
    retval = my_rise_trans(JDNDaysUT, -1, ObjectName, eventtype, helflag, dgeo, datm, tret, serr);
  return retval;
}

/*###################################################################
' JDNDaysUT [Days]
' actual [0= approximation, 1=actual]
' SunRA [deg]
*/
static double SunRA(double JDNDaysUT, int32 helflag, char *serr)
{
  int imon, iday, iyar, calflag = SE_GREG_CAL;
  double dut;
  static TLS double tjdlast;
  static TLS double ralast;
  helflag += 0; /* statement prevents compiler warning */
  *serr = '\0';
  if (JDNDaysUT == tjdlast)
    return ralast;
#ifndef SIMULATE_VICTORVB
  if (1) { /*helflag & SE_HELFLAG_HIGH_PRECISION) {*/
    double tjd_tt;
    double x[6];
    int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
    int32 iflag = epheflag | SEFLG_EQUATORIAL;
    iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
    tjd_tt = JDNDaysUT + swe_deltat_ex(JDNDaysUT, epheflag, serr);
    if (swe_calc(tjd_tt, SE_SUN, iflag, x, serr) != ERR) {
      ralast = x[0];
      tjdlast = JDNDaysUT;
      return ralast;
    }
  }
#endif
  swe_revjul(JDNDaysUT, calflag, &iyar, &imon, &iday, &dut); /* this seems to be much faster than calling swe_revjul() ! Note: only because SunRA is called 1000s of times */
  tjdlast = JDNDaysUT;
  ralast = swe_degnorm((imon + (iday - 1) / 30.4 - 3.69) * 30);
  /*ralast = (DatefromJDut(JDNDaysUT, 2) + (DatefromJDut(JDNDaysUT, 3) - 1) / 30.4 - 3.69) * 30;*/
  return ralast;
}

/*###################################################################
' Temp [C]
' Kelvin [K]
*/
static double Kelvin(double Temp)
{
  /*' http://en.wikipedia.org/wiki/Kelvin*/
  return Temp + C2K;
}

/*###################################################################
' AppAlt [deg]
' TempE [C]
' PresE [mbar]
' TopoAltitudefromAppAlt [deg]
*/
static double TopoAltfromAppAlt(double AppAlt, double TempE, double PresE)
{
  double R = 0;
  double retalt = 0;
  if (AppAlt >= LowestAppAlt) {
    if (AppAlt > 17.904104638432)
      R = 0.97 / tan(AppAlt * DEGTORAD);
    else
      R = (34.46 + 4.23 * AppAlt + 0.004 * AppAlt * AppAlt) / (1 + 0.505 * AppAlt + 0.0845 * AppAlt * AppAlt);
    R = (PresE - 80) / 930 / (1 + 0.00008 * (R + 39) * (TempE - 10)) * R;
    retalt = AppAlt - R * Min2Deg;
  } else {
    retalt = AppAlt;
  }
  return retalt;
}

/*###################################################################
' TopoAlt [deg]
' TempE [C]
' PresE [mbar]
' AppAltfromTopoAlt [deg]
' call this instead of swe_azalt(), because it is faster (lower precision
' is required)
*/
static double AppAltfromTopoAlt(double TopoAlt, double TempE, double PresE, int32 helflag)
{
  /* using methodology of Newtown derivatives (analogue to what Swiss Emphemeris uses)*/
  int i, nloop = 2;
  double newAppAlt = TopoAlt;
  double newTopoAlt = 0.0;
  double oudAppAlt = newAppAlt;
  double oudTopoAlt = newTopoAlt;
  double verschil, retalt;
  if (helflag & SE_HELFLAG_HIGH_PRECISION)
    nloop = 5;
  for (i = 0; i <= nloop; i++) {
    newTopoAlt = newAppAlt - TopoAltfromAppAlt(newAppAlt, TempE, PresE);
    /*newTopoAlt = newAppAlt - swe_refrac(newAppAlt, PresE, TempE, SE_CALC_APP_TO_TRUE);*/
    verschil = newAppAlt - oudAppAlt;
    oudAppAlt = newTopoAlt - oudTopoAlt - verschil;
    if ((verschil != 0) && (oudAppAlt != 0))
      verschil = newAppAlt - verschil * (TopoAlt + newTopoAlt - newAppAlt) / oudAppAlt;
    else
      verschil = TopoAlt + newTopoAlt;
    oudAppAlt = newAppAlt;
    oudTopoAlt = newTopoAlt;
    newAppAlt = verschil;
  }
  retalt = TopoAlt + newTopoAlt;
  if (retalt < LowestAppAlt)
    retalt = TopoAlt;
  return retalt;
}

/*###################################################################
' TopoAlt [deg]
' TopoDecl [deg]
' Lat [deg]
' HourAngle [hour]
*/
static double HourAngle(double TopoAlt, double TopoDecl, double Lat)
{
  double Alti = TopoAlt * DEGTORAD;
  double decli = TopoDecl * DEGTORAD;
  double Lati = Lat * DEGTORAD;
  double ha = (sin(Alti) - sin(Lati) * sin(decli)) / cos(Lati) / cos(decli);
  if (ha < -1) ha = -1; 
  if (ha > 1) ha = 1;
  /* from http://star-www.st-and.ac.uk/~fv/webnotes/chapt12.htm*/
  return acos(ha) / DEGTORAD / 15.0;
}

/*###################################################################
' JDNDaysUT [Days]
' dgeo [array: longitude, latitude, eye height above sea m]
' TempE [C]
' PresE [mbar]
' ObjectName [-]
' Angle (0 = TopoAlt, 1 = Azi, 2=Topo Declination, 3=Topo Rectascension, 4=AppAlt,5=Geo Declination, 6=Geo Rectascension)
' ObjectLoc [deg]
 */
static int32 ObjectLoc(double JDNDaysUT, double *dgeo, double *datm, char *ObjectName, int32 Angle, int32 helflag, double *dret, char *serr)
{
  double x[6], xin[3], xaz[3], tjd_tt;
  int32 Planet;
  int32 epheflag;
  int32 iflag = SEFLG_EQUATORIAL;
  epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  iflag |= epheflag;
  if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
    iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
  if (Angle < 5) iflag = iflag | SEFLG_TOPOCTR;
  if (Angle == 7) Angle = 0;
  tjd_tt = JDNDaysUT + swe_deltat_ex(JDNDaysUT, epheflag, serr);
  Planet = DeterObject(ObjectName);
  if (Planet != -1) {
    if (swe_calc(tjd_tt, Planet, iflag, x, serr) == ERR)
      return ERR;
  } else {
    if (call_swe_fixstar(ObjectName, tjd_tt, iflag, x, serr) == ERR)
      return ERR;
  }
  if (Angle == 2 ||  Angle == 5) {
    *dret = x[1];
  } else {
    if (Angle == 3 || Angle == 6) {
      *dret = x[0];
    } else {
      xin[0] = x[0];
      xin[1] = x[1];
      swe_azalt(JDNDaysUT, SE_EQU2HOR, dgeo, datm[0], datm[1], xin, xaz);
      if (Angle == 0)
	*dret = xaz[1];
      if (Angle == 4)
	*dret = AppAltfromTopoAlt(xaz[1], datm[0], datm[1], helflag);
      if (Angle == 1) {
	xaz[0] += 180;
	if (xaz[0] >= 360)
	  xaz[0] -= 360;
        *dret = xaz[0];
      }
    }
  }
  return OK;
}

/*###################################################################
' JDNDaysUT [Days]
' dgeo [array: longitude, latitude, eye height above sea m]
' TempE [C]
' PresE [mbar]
' ObjectName [-]
' Angle (0 = TopoAlt, 1 = Azi, 2=Topo Declination, 3=Topo Rectascension, 4=AppAlt,5=Geo Declination, 6=Geo Rectascension)
' ObjectLoc [deg]
 */
static int32 azalt_cart(double JDNDaysUT, double *dgeo, double *datm, char *ObjectName, int32 helflag, double *dret, char *serr)
{
  double x[6], xin[3], xaz[3], tjd_tt;
  int32 Planet;
  int32 epheflag;
  int32 iflag = SEFLG_EQUATORIAL;
  epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  iflag |= epheflag;
  if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
    iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
  iflag = iflag | SEFLG_TOPOCTR;
  tjd_tt = JDNDaysUT + swe_deltat_ex(JDNDaysUT, epheflag, serr);
  Planet = DeterObject(ObjectName);
  if (Planet != -1) {
    if (swe_calc(tjd_tt, Planet, iflag, x, serr) == ERR)
      return ERR;
  } else {
    if (call_swe_fixstar(ObjectName, tjd_tt, iflag, x, serr) == ERR)
      return ERR;
  }
  xin[0] = x[0];
  xin[1] = x[1];
  swe_azalt(JDNDaysUT, SE_EQU2HOR, dgeo, datm[0], datm[1], xin, xaz);
  dret[0] = xaz[0];
  dret[1] = xaz[1]; /* true altitude */
  dret[2] = xaz[2]; /* apparent altitude */
  /* also return cartesian coordinates, for apparent altitude */
  xaz[1] = xaz[2];
  xaz[2] = 1;
  swi_polcart(xaz, xaz);
  dret[3] = xaz[0];
  dret[4] = xaz[1];
  dret[5] = xaz[2];
  return OK;
}

/*###################################################################
' LatA [rad]
' LongA [rad]
' LatB [rad]
' LongB [rad]
' DistanceAngle [rad]
*/
static double DistanceAngle(double LatA, double LongA, double LatB, double LongB)
{
  double dlon = LongB - LongA;
  double dlat = LatB - LatA;
  /* Haversine formula
   * http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html
   * R.W. Sinnott, Virtues of the Haversine, Sky and Telescope, vol. 68, no. 2, 1984, p. 159
   */
  double sindlat2 = sin(dlat / 2);
  double sindlon2 = sin(dlon / 2);
  double corde = sindlat2 * sindlat2 + cos(LatA) * cos(LatB) * sindlon2 *sindlon2;
  if (corde > 1) corde = 1;
  return 2 * asin(sqrt(corde));
}

/*###################################################################
' heighteye [m]
' TempS [C]
' RH [%]
' kW [-]
*/
static double kW(double HeightEye, double TempS, double RH)
{
  /* From Schaefer , Archaeoastronomy, XV, 2000, page 128*/
  double WT = 0.031;
  WT *= 0.94 * (RH / 100.0) * exp(TempS / 15) * exp(-1 * HeightEye / scaleHwater);
  return WT;
}

/*###################################################################
' JDNDaysUT [-]
' AltS [deg]
' lat [deg]
' kOZ [-]
*/
static double kOZ(double AltS, double sunra, double Lat)
{
  double CHANGEKO, OZ, LT, kOZret;
  static TLS double koz_last, alts_last, sunra_last;
  double altslim = 0;
  if (AltS == alts_last && sunra == sunra_last)
    return koz_last;
  alts_last = AltS; sunra_last = sunra;
  OZ = 0.031;
  LT = Lat * DEGTORAD;
  /* From Schaefer , Archaeoastronomy, XV, 2000, page 128*/
  kOZret = OZ * (3.0 + 0.4 * (LT * cos(sunra * DEGTORAD) - cos(3 * LT))) / 3.0;
  /* depending on day/night vision (altitude of sun < start astronomical twilight), KO changes from 100% to 30%
   * see extinction section of Vistas in Astronomy page 343*/
  altslim = -AltS - 12;
  if (altslim < 0)
    altslim = 0;
  CHANGEKO = (100 - 11.6 * mymin(6, altslim)) / 100;
if ((0)) {
  static int a = 0;
  if (a == 0)
    printf("bsk=%f %f\n", kOZret, AltS);
  a = 1;
}
  koz_last = kOZret * CHANGEKO;
  return koz_last;
}

/*###################################################################
' AltS [deg]
' heighteye [m]
' kR [-]
*/
static double kR(double AltS, double HeightEye)
{
  /* depending on day/night vision (altitude of sun < start astronomical twilight),
   * lambda eye sensibility changes
   * see extinction section of Vistas in Astronomy page 343*/
  double CHANGEK, LAMBDA;
  double val = -AltS - 12;
  if (val < 0) val = 0;
  if (val > 6) val = 6;
  /*CHANGEK = (1 - 0.166667 * Min(6, Max(-AltS - 12, 0)));*/
  CHANGEK = (1 - 0.166667 * val );
  LAMBDA = 0.55 + (CHANGEK - 1) * 0.04;
  /* From Schaefer , Archaeoastronomy, XV, 2000, page 128 */
  return 0.1066 * exp(-1 * HeightEye / scaleHrayleigh) * pow(LAMBDA / 0.55 , -4);
}

static int Sgn(double x)
{
  if (x < 0) 
    return -1;
  return 1;
}

/*###################################################################
' JDNDaysUT [-]
' AltS [deg]
' lat [deg]
' heighteye [m]
' TempS [C]
' RH [%]
' VR [km]
' ka [-]
*/
static double ka(double AltS, double sunra, double Lat, double HeightEye, double TempS, double RH, double VR, char *serr)
{
  double CHANGEKA, LAMBDA, BetaVr, Betaa, kaact;
  double SL = Sgn(Lat);
  /* depending on day/night vision (altitude of sun < start astronomical twilight),
   * lambda eye sensibility changes
   * see extinction section of Vistas in Astronomy page 343 */
  static TLS double alts_last, sunra_last, ka_last;
  if (AltS == alts_last && sunra == sunra_last)
    return ka_last;
  alts_last = AltS; sunra_last = sunra;
  CHANGEKA = (1 - 0.166667 * mymin(6, mymax(-AltS - 12, 0)));
  LAMBDA = 0.55 + (CHANGEKA - 1) * 0.04;
  if (VR != 0) {
    if (VR >= 1) {
      /* Visbility range from http://www1.cs.columbia.edu/CAVE/publications/pdfs/Narasimhan_CVPR03.pdf
       * http://www.icao.int/anb/SG/AMOSSG/meetings/amossg3/wp/SN11Rev.pdf where MOR=2.995/ke
       * factor 1.3 is the relation between "prevailing visibility" and 
       * meteorological range was derived by Koshmeider in the 1920's */
      BetaVr = 3.912 / VR;
      Betaa = BetaVr - (kW(HeightEye, TempS, RH) / scaleHwater + kR(AltS, HeightEye) / scaleHrayleigh) * 1000 * astr2tau;
      kaact = Betaa * scaleHaerosol / 1000 * tau2astr;
      if (kaact < 0) {
	if (serr != NULL)
	  strcpy(serr, "The provided Meteorological range is too long, when taking into acount other atmospheric parameters"); /* is a warning */
        /* return 0; * return "#HIGHVR"; */
      }
    } else {
      kaact = VR - kW(HeightEye, TempS, RH) - kR(AltS, HeightEye) - kOZ(AltS, sunra, Lat);
      if (kaact < 0) {
	if (serr != NULL)
	  strcpy(serr, "The provided atmosphic coeefficent (ktot) is too low, when taking into acount other atmospheric parameters"); /* is a warning */
        /* return 0; * "#LOWktot"; */
      }
    }
  } else {
    /* From Schaefer , Archaeoastronomy, XV, 2000, page 128 */
#ifdef SIMULATE_VICTORVB
    if (RH <= 0.00000001) RH = 0.00000001;
    if (RH >= 99.99999999) RH = 99.99999999;
#endif
    kaact = 0.1 * exp(-1 * HeightEye / scaleHaerosol) * pow(1 - 0.32 / log(RH / 100.0), 1.33) * (1 + 0.33 * SL * sin(sunra * DEGTORAD));
    kaact = kaact * pow(LAMBDA / 0.55, -1.3);
  }
  ka_last = kaact;
  return kaact;
}

/*###################################################################
' JDNDaysUT [-]
' AltS [deg]
' lat [deg]
' heighteye [m]
' TempS [C]
' RH [%]
' VR [km]
' ExtType [0=ka,1=kW,2=kR,3=kOZ,4=ktot]
' kt [-]
*/
static double kt(double AltS, double sunra, double Lat, double HeightEye, double TempS, double RH, double VR, int32 ExtType, char *serr)
{
  double kRact = 0;
  double kWact = 0;
  double kOZact = 0;
  double kaact = 0;
  if (ExtType == 2 || ExtType == 4)
    kRact = kR(AltS, HeightEye);
  if (ExtType == 1 || ExtType == 4)
    kWact = kW(HeightEye, TempS, RH);
  if (ExtType == 3 || ExtType == 4)
    kOZact = kOZ(AltS, sunra, Lat);
  if (ExtType == 0 || ExtType == 4)
    kaact = ka(AltS, sunra, Lat, HeightEye, TempS, RH, VR, serr);
  if (kaact < 0)
    kaact = 0;
  return kWact + kRact + kOZact + kaact;
}

/*###################################################################
' AppAlt0 [deg]
' PresS [mbar]
' Airmass [??]
*/
static double Airmass(double AppAltO, double Press)
{
  double airm, zend;
  zend = (90 - AppAltO) * DEGTORAD;
  if (zend > PI / 2)
    zend = PI / 2;
  airm = 1 / (cos(zend) + 0.025 * exp(-11 * cos(zend)));
  return Press / 1013 * airm;
}

/*###################################################################
' scaleH '[m]
' zend [rad]
' PresS [mbar]
' Xext [-]
*/
static double Xext(double scaleH, double zend, double Press)
{
  return Press / 1013.0 / (cos(zend) + 0.01 * sqrt(scaleH / 1000.0) * exp(-30.0 / sqrt(scaleH / 1000.0) * cos(zend)));
}

/*###################################################################
' scaleH '[m]
' zend [rad]
' PresS [mbar]
' Xlay [-]
*/
static double Xlay(double scaleH, double zend, double Press)
{
  /*return Press / 1013.0 / sqrt(1.0 - pow(sin(zend) / (1.0 + (scaleH / Ra)), 2));*/
  double a = sin(zend) / (1.0 + (scaleH / Ra));
  return Press / 1013.0 / sqrt(1.0 - a * a);
}

/*###################################################################
' Meteorological formula
'###################################################################
' TempS [C]
' HeightEye [m]
' TempEfromTempS [C]
*/
static double TempEfromTempS(double TempS, double HeightEye, double Lapse)
{
  return TempS - Lapse * HeightEye;
}

/*###################################################################
' TempS [C]
' PresS [mbar]
' HeightEye [m]
' PresEfromPresS [mbar]
*/
static double PresEfromPresS(double TempS, double Press, double HeightEye)
{
  return Press * exp(-9.80665 * 0.0289644 / (Kelvin(TempS) + 3.25 * HeightEye / 1000) / 8.31441 * HeightEye);
}

/*###################################################################
' AltO [deg]
' JDNDaysUT [-]
' AltS [deg]
' lat [deg]
' heighteye [m]
' TempS [C]
' PresS [mbar]
' RH [%]
' VR [km]
' Deltam [-]
*/
static double Deltam(double AltO, double AltS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
{
  double zend, xR, XW, Xa, XOZ;
  double PresE = PresEfromPresS(datm[1], datm[0], HeightEye);
  double TempE = TempEfromTempS(datm[1], HeightEye, LapseSA);
  double AppAltO = AppAltfromTopoAlt(AltO, TempE, PresE, helflag);
  double deltam;
  static TLS double alts_last, alto_last, sunra_last, deltam_last;
  if (AltS == alts_last && AltO == alto_last && sunra == sunra_last)
    return deltam_last;
  alts_last = AltS; alto_last = AltO; sunra_last = sunra;
  if (staticAirmass == 0) {
    zend = (90 - AppAltO) * DEGTORAD;
    if (zend > PI / 2)
      zend = PI / 2;
    /* From Schaefer , Archaeoastronomy, XV, 2000, page 128*/
    xR = Xext(scaleHrayleigh, zend, datm[0]);
    XW = Xext(scaleHwater, zend, datm[0]);
    Xa = Xext(scaleHaerosol, zend, datm[0]);
    XOZ = Xlay(scaleHozone, zend, datm[0]);
    deltam = kR(AltS, HeightEye) * xR + kt(AltS, sunra, Lat, HeightEye, datm[1], datm[2], datm[3], 0, serr) * Xa + kOZ(AltS, sunra, Lat) * XOZ + kW(HeightEye, datm[1], datm[2]) * XW;
  } else {
    deltam = kt(AltS, sunra, Lat, HeightEye, datm[1], datm[2], datm[3], 4, serr) * Airmass(AppAltO, datm[0]);
  }
  deltam_last = deltam;
  return deltam;
}

/*###################################################################
' AltO [deg]
' JDNDaysUT [-]
' AltS [deg]
' lat [deg]
' heighteye [m]
' TempS [C]
' PresS [mbar]
' RH [%]
' VR [km]
' Bn [nL]
*/
static double Bn(double AltO, double JDNDayUT, double AltS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
{
  double PresE = PresEfromPresS(datm[1], datm[0], HeightEye);
  double TempE = TempEfromTempS(datm[1], HeightEye, LapseSA);
  double AppAltO = AppAltfromTopoAlt(AltO, TempE, PresE, helflag);
  double zend, YearB, MonthB, DayB, Bna, kX, Bnb;
  double B0 = 0.0000000000001, dut;
  int iyar, imon, iday;
  /* Below altitude of 10 degrees, the Bn stays the same (see page 343 Vistas in Astronomy) */
  if (AppAltO < 10)
    AppAltO = 10;
  zend = (90 - AppAltO) * DEGTORAD;
  /* From Schaefer , Archaeoastronomy, XV, 2000, page 128 and adjusted for sunspot period*/
  /*YearB = DatefromJDut(JDNDayUT, 1);
    MonthB = DatefromJDut(JDNDayUT, 2);
    DayB = DatefromJDut(JDNDayUT, 3);*/
  swe_revjul(JDNDayUT, SE_GREG_CAL, &iyar, &imon, &iday, &dut); 
  YearB = iyar; MonthB = imon; DayB = iday;
  Bna = B0 * (1 + 0.3 * cos(6.283 * (YearB + ((DayB - 1) / 30.4 + MonthB - 1) / 12 - 1990.33) / 11.1));
  kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
  /* From Schaefer , Archaeoastronomy, XV, 2000, page 129 */
  Bnb = Bna * (0.4 + 0.6 / sqrt(1 - 0.96 * pow(sin(zend), 2))) * pow(10, -0.4 * kX);
  return mymax(Bnb, 0) * erg2nL;
}

/*###################################################################
' JDNDaysUT [-]
' dgeo [array: longitude, latitude, eye height above sea m]
' TempE [C]
' PresE [mbar]
' ObjectName [-]
' Magnitude [-]
*/
static int32 Magnitude(double JDNDaysUT, double *dgeo, char *ObjectName, int32 helflag, double *dmag, char *serr)
{
  double x[20];
  int32 Planet, iflag, epheflag;
  epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  *dmag = -99.0;
  Planet = DeterObject(ObjectName);
  iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
  if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
    iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
  if (Planet != -1) {
    /**dmag = Phenomena(JDNDaysUT, Lat, Longitude, HeightEye, TempE, PresE, ObjectName, 4);*/
    swe_set_topo(dgeo[0], dgeo[1], dgeo[2]);
    if (swe_pheno_ut(JDNDaysUT, Planet, iflag, x, serr) == ERR)
      return ERR;
    *dmag = x[4];
  } else {
    if (call_swe_fixstar_mag(ObjectName, dmag, serr) == ERR)
      return ERR;
  }
  return OK;
}

#if 0
static int32 fast_magnitude(double tjd, double *dgeo, char *ObjectName, int32 helflag, double *dmag, char *serr)
{
  int32 retval = OK, ipl, ipli;
  double dtjd;
  static TLS double tjdsv[3];
  static TLS double dmagsv[3];
  static TLS int32 helflagsv[3];
  ipl = DeterObject(ObjectName);
  ipli = ipl;
  if (ipli > SE_MOON) 
    ipli = 2;
  dtjd = tjd - tjdsv[ipli];
  if (tjdsv[ipli] != 0 && helflag == helflagsv[ipli] && fabs(dtjd) < 5.0 / 1440.0) {
    *dmag = dmagsv[ipli];
  } else {
    retval = Magnitude(tjd, dgeo, ObjectName, helflag, dmag, serr);
    tjdsv[ipli] = tjd;
    helflagsv[ipli] = helflag;
    dmagsv[ipli] = *dmag;
  }
  return retval;
}
#endif

/*###################################################################
' dist [km]
' phasemoon [-]
' MoonsBrightness [-]
*/
static double MoonsBrightness(double dist, double phasemoon)
{
  double log10 = 2.302585092994;
  /*Moon's brightness changes with distance: http://hem.passagen.se/pausch/comp/ppcomp.html#15 */
  return -21.62 + 5 * log(dist / (Ra / 1000)) / log10 + 0.026 * fabs(phasemoon) + 0.000000004 * pow(phasemoon, 4);
}

/*###################################################################
' AltM [deg]
' AziM [deg]
' AziS [deg]
' MoonPhase [deg]
*/
static double MoonPhase(double AltM, double AziM, double AltS, double AziS)
{
  double AltMi = AltM * DEGTORAD;
  double AltSi = AltS * DEGTORAD;
  double AziMi = AziM * DEGTORAD;
  double AziSi = AziS * DEGTORAD;
  double MoonAvgPar = 0.95;
// return 180 - acos(cos(AziSi - AziMi) * cos(AltMi + MoonAvgPar * DEGTORAD) * cos(AltSi) + sin(AltSi) * sin(AltMi + MoonAvgPar * DEGTORAD)) / DEGTORAD;
  return 180 - acos(cos(AziSi - AziMi - MoonAvgPar * DEGTORAD) * cos(AltMi + MoonAvgPar * DEGTORAD) * cos(AltSi) + sin(AltSi) * sin(AltMi + MoonAvgPar * DEGTORAD)) / DEGTORAD;
}

/*###################################################################
' Pressure [mbar]
*/
static double Bm(double AltO, double AziO, double AltM, double AziM, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
{
  double M0 = -11.05;
  double Bm = 0;
  double RM, kXM, kX, C3, FM, phasemoon, MM;
  double lunar_radius = 0.25 * DEGTORAD;
  AS_BOOL object_is_moon = FALSE;
  if (AltO == AltM && AziO == AziM)
    object_is_moon = TRUE;
  if (AltM > -0.26 && !object_is_moon) { // second condition added by Dieter, SE2.06
  /* moon only adds light when (partly) above horizon
   * From Schaefer , Archaeoastronomy, XV, 2000, page 129*/
    RM = DistanceAngle(AltO * DEGTORAD, AziO * DEGTORAD, AltM * DEGTORAD, AziM * DEGTORAD) / DEGTORAD;
    if (RM <= lunar_radius) // addition by Dieter for objects behind the Moon, SE2.06
      RM = lunar_radius;
    kXM = Deltam(AltM, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
    kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
    C3 = pow(10, -0.4 * kXM);
    FM = (62000000.0) / RM / RM + pow(10, 6.15 - RM / 40) + pow(10, 5.36) * (1.06 + pow(cos(RM * DEGTORAD), 2));
    Bm = FM * C3 + 440000 * (1 - C3);
    phasemoon = MoonPhase(AltM, AziM, AltS, AziS);
    MM = MoonsBrightness(MoonDistance, phasemoon);
    Bm = Bm * pow(10, -0.4 * (MM - M0 + 43.27));
    Bm = Bm * (1 - pow(10, -0.4 * kX));
  }
  Bm = mymax(Bm, 0) * erg2nL;
  return Bm;
}

/*###################################################################
' Pressure [mbar]
*/
static double Btwi(double AltO, double AziO, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
{
double M0 = -11.05;
double MS = -26.74;
double PresE = PresEfromPresS(datm[1], datm[0], HeightEye);
double TempE = TempEfromTempS(datm[1], HeightEye, LapseSA);
double AppAltO = AppAltfromTopoAlt(AltO, TempE, PresE, helflag);
double ZendO = 90 - AppAltO;
double RS = DistanceAngle(AltO * DEGTORAD, AziO * DEGTORAD, AltS * DEGTORAD, AziS * DEGTORAD) / DEGTORAD;
double kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
double k = kt(AltS, sunra, Lat, HeightEye, datm[1], datm[2], datm[3], 4, serr);
/* From Schaefer , Archaeoastronomy, XV, 2000, page 129*/
double Btwi = pow(10, -0.4 * (MS - M0 + 32.5 - AltS - (ZendO / (360 * k))));
Btwi = Btwi * (100 / RS) * (1 - pow(10, -0.4 * kX));
Btwi = mymax(Btwi, 0) * erg2nL;
return Btwi;
}

/*###################################################################
' Pressure [mbar]
2300 REM  Daylight brightness
2310 C4=10.0^(-.4*K(I)*XS)
2320 FS=6.2E+07*(RS^-2)+(10^(6.15-RS/40))
2330 FS=FS+(10^5.36)*(1.06+((COS(RS*RD))^2))
2340 BD=10^(-.4*(MS(I)-MO(I)+43.27))
2350 BD=BD*(1-10^(-.4*K(I)*X))
2360 BD=BD*(FS*C4+440000.0*(1-C4))
*/
static double Bday(double AltO, double AziO, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
{
  double M0 = -11.05;
  double MS = -26.74;
  double RS = DistanceAngle(AltO * DEGTORAD, AziO * DEGTORAD, AltS * DEGTORAD, AziS * DEGTORAD) / DEGTORAD;
  double kXS = Deltam(AltS, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
  double kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
  /* From Schaefer , Archaeoastronomy, XV, 2000, page 129*/
  double C4 = pow(10, -0.4 * kXS);
  double FS = (62000000.0) / RS / RS + pow(10, (6.15 - RS / 40)) + pow(10, 5.36) * (1.06 + pow(cos(RS * DEGTORAD), 2));
  double Bday = FS * C4 + 440000.0 * (1 - C4);
  Bday = Bday * pow(10, (-0.4 * (MS - M0 + 43.27)));
  Bday = Bday * (1 - pow(10, -0.4 * kX));
  Bday = mymax(Bday, 0) * erg2nL;
  return Bday;
}

/*###################################################################
' Value [nL]
' PresS [mbar]
' Bcity [nL]
*/
static double Bcity(double Value, double Press)
{
  double Bcity = Value;
  Press += 0.0; /* unused; statement prevents compiler warning */
  Bcity = mymax(Bcity, 0);
  return Bcity;
}

/*###################################################################
' Pressure [mbar]
*/
static double Bsky(double AltO, double AziO, double AltM, double AziM, double JDNDaysUT, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
{
  double Bsky = 0;
  if (AltS < -3) {
    Bsky += Btwi(AltO, AziO, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr);
  } else {
    if (AltS > 4) {
      Bsky += Bday(AltO, AziO, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr);
    } else {
      Bsky += mymin(Bday(AltO, AziO, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr), Btwi(AltO, AziO, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr));
if ((0)) {
  static int a = 0;
  if (a == 0)
    printf("bsk=%f\n", Bsky);
  a = 1;
}
    }
  }
  /* if max. Bm [1E7] <5% of Bsky don't add Bm*/
  if (Bsky < 200000000.0) 
    Bsky += Bm(AltO, AziO, AltM, AziM, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr);
  if (AltS <= 0)
    Bsky += Bcity(0, datm[0]);
  /* if max. Bn [250] <5% of Bsky don't add Bn*/
  if (Bsky < 5000)
    Bsky = Bsky + Bn(AltO, JDNDaysUT, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
  /* if max. Bm [1E7] <5% of Bsky don't add Bm*/
  return Bsky;
}

/* default handling:
 * 1. datm (atmospheric conditions):
 * datm consists of 
 *     [0]  atmospheric pressure
 *     [1]  temperature
 *     [2]  relative humidity
 *     [3]  extinction coefficient
 * In order to get default values for [0..2], set datm[0] = 0.
 * Default values for [1-2] are only provided if [0] == 0.
 * [3] defaults outside this function, depending on [0-2].
 * 
 * 2. dobs (observer definition):
 *     [0]  age (default 36)
 *     [1]  Snellen ratio or visual acuity of observer (default 1)
 */
static void default_heliacal_parameters(double *datm, double *dgeo, double *dobs, int helflag)
{
  int i;
  if (datm[0] <= 0) {
    /* estimate atmospheric pressure, according to the
     * International Standard Atmosphere (ISA) */
    datm[0] = 1013.25 * pow(1 - 0.0065 * dgeo[2] / 288, 5.255);
    /* temperature */
    if (datm[1] == 0)
      datm[1] = 15 - 0.0065 * dgeo[2];
    /* relative humidity, independent of atmospheric pressure and altitude */
    if (datm[2] == 0)
      datm[2] = 40;
    /* note: datm[3] / VR defaults outside this function */
  } else {
#ifndef SIMULATE_VICTORVB
    if (datm[2] <= 0.00000001) datm[2] = 0.00000001;
    if (datm[2] >= 99.99999999) datm[2] = 99.99999999;
#endif
  }
  /* age of observer */
  if (dobs[0] == 0)
    dobs[0] = 36;
  /* SN Snellen factor of the visual acuity of the observer */
  if (dobs[1] == 0)
    dobs[1] = 1;
  if (!(helflag & SE_HELFLAG_OPTICAL_PARAMS)) {
    for (i = 2; i <= 5; i++)
      dobs[i] = 0;
  }
  /* OpticMagn undefined -> use eye */
  if (dobs[3] == 0) {
    dobs[2] = 1; /* Binocular = 1 */
    dobs[3] = 1; /* OpticMagn = 1: use eye */
    /* dobs[4] and dobs[5] (OpticDia and OpticTrans) will be defaulted in 
     * OpticFactor() */
  }
}

/*###################################################################
' age [Year]
' SN [-]
' AltO [deg]
' AziO [deg]
' AltM [deg]
' AziM [deg]
' MoonDistance [km]
' JDNDaysUT [-]
' AltS [deg]
' AziS [deg]
' lat [deg]
' heighteye [m]
' TempS [C]
' PresS [mbar]
' RH [%]
' VR [km]
' VisLimMagn [-]
*/
static double VisLimMagn(double *dobs, double AltO, double AziO, double AltM, double AziM, double JDNDaysUT, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, int32 *scotopic_flag, char *serr)
{
  double C1, C2, Th, kX, Bsk, CorrFactor1, CorrFactor2;
  double log10 = 2.302585092994;
  AS_BOOL is_scotopic = FALSE;
  /*double Age = dobs[0];*/
  /*double SN = dobs[1];*/
  Bsk = Bsky(AltO, AziO, AltM, AziM, JDNDaysUT, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr);
  /* Schaefer, Astronomy and the limits of vision, Archaeoastronomy, 1993 Verder:*/
  kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
if ((0)) {
  static int a = 0;
  if (a == 0)
    printf("bsk=%f, kx=%f\n", Bsk, kX);
  a = 1;
}
  /* influence of age*/
  /*Fa = mymax(1, pow(p(23, Bsk) / p(Age, Bsk), 2)); */
  CorrFactor1 = OpticFactor(Bsk, kX, dobs, JDNDaysUT, "", 1, helflag);
  CorrFactor2 = OpticFactor(Bsk, kX, dobs, JDNDaysUT, "", 0, helflag);
  //if (Bsk < BNIGHT)
  if (Bsk < 1645) /* use this function for BNIGHT to make the function continuous */
    is_scotopic = TRUE;
  if (helflag & SE_HELFLAG_VISLIM_PHOTOPIC) 
    is_scotopic = FALSE;
  if (helflag & SE_HELFLAG_VISLIM_SCOTOPIC) 
    is_scotopic = TRUE;
  /* From Schaefer , Archaeoastronomy, XV, 2000, page 129*/
  if (is_scotopic) {
    C1 = 1.5848931924611e-10; /*pow(10, -9.8);*/ /* C1 = 10 ^ (-9.8);*/
    C2 = 0.012589254117942; /*pow(10, -1.9);*/ /* C2 = 10 ^ (-1.9);*/
    if (scotopic_flag != NULL) 
      *scotopic_flag = 1;
  } else {
    C1 = 4.4668359215096e-9; /*pow(10, -8.35);*/ /* C1 = 10 ^ (-8.35);*/
    C2 = 1.2589254117942e-6; /*pow(10, -5.9);*/ /* C2 = 10 ^ (-5.9);*/
    if (scotopic_flag != NULL) 
      *scotopic_flag = 0;
  }
  if (scotopic_flag != NULL) {
    if (BNIGHT * BNIGHT_FACTOR > Bsk && BNIGHT / BNIGHT_FACTOR < Bsk)
      *scotopic_flag |= 2;
  }
  /*Th = C1 * pow(1 + sqrt(C2 * Bsk), 2) * Fa;*/
  /*Bsk = Bsk / CorrFactor1;*/
  Bsk = Bsk * CorrFactor1;
  Th = C1 * pow(1 + sqrt(C2 * Bsk), 2) * CorrFactor2;
#if SWEHEL_DEBUG
  fprintf(stderr, "Bsk=%f, ", Bsk);
  fprintf(stderr, "kX =%f, ", kX);
  fprintf(stderr, "Th =%f, ", Th);
  fprintf(stderr, "CorrFactor1=%f, ", CorrFactor1);
  fprintf(stderr, "CorrFactor2=%f\n", CorrFactor2);
#endif
  /* Visual limiting magnitude of point source*/
#if 0
  if (SN <= 0.00000001)
    SN = 0.00000001;
  return -16.57 - 2.5 * (log(Th) / log10) - kX + 5.0 * (log(SN) / log10);*/
#endif
  return -16.57 - 2.5 * (log(Th) / log10);
}

/* tolower star name, but not Bayer designation */
static char *tolower_string_star(char *str)
{
  char *sp;
  for (sp = str; *sp != '\0' && *sp != ','; sp++)
    *sp = tolower(*sp);
  return str;
}

/* Limiting magnitude in dark skies 
 * for information about input parameters, see function swe_heliacal_ut().
 *
 * function returns:
 * -1   Error
 * -2   Object is below horizon
 *  0   OK, photopic vision
 *  |1  OK, scotopic vision
 *  |2  OK, near limit photopic/scotopic
*/
int32 CALL_CONV swe_vis_limit_mag(double tjdut, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, double *dret, char *serr)
{
  int32 retval = OK, i, scotopic_flag = 0;
  double AltO, AziO, AltM, AziM, AltS, AziS;
  double sunra;
  for (i = 0; i < 7; i++)
    dret[i] = 0;
  tolower_string_star(ObjectName);
  if (DeterObject(ObjectName) == SE_SUN) {
    if (serr != NULL) {
      strcpy(serr, "it makes no sense to call swe_vis_limit_mag() for the Sun");
    }
    return ERR;
  }
  swi_set_tid_acc(tjdut, helflag, 0, serr);
  sunra = SunRA(tjdut, helflag, serr);
  default_heliacal_parameters(datm, dgeo, dobs, helflag);
  swe_set_topo(dgeo[0], dgeo[1], dgeo[2]);
  if (ObjectLoc(tjdut, dgeo, datm, ObjectName, 0, helflag, &AltO, serr) == ERR)
    return ERR;
  if (AltO < 0) {
    if (serr != NULL)
      strcpy(serr, "object is below local horizon");
    *dret = -100;
    return -2;
  }
  if (ObjectLoc(tjdut, dgeo, datm, ObjectName, 1, helflag, &AziO, serr) == ERR)
    return ERR;
  if (helflag & SE_HELFLAG_VISLIM_DARK) {
    AltS = -90;
    AziS = 0;
  } else {
    if (ObjectLoc(tjdut, dgeo, datm, "sun", 0, helflag, &AltS, serr) == ERR)
      return ERR;
    if (ObjectLoc(tjdut, dgeo, datm, "sun", 1, helflag, &AziS, serr) == ERR)
      return ERR;
  }
  if (strncmp(ObjectName, "moon", 4) == 0 ||
      (helflag & SE_HELFLAG_VISLIM_DARK) ||
      (helflag & SE_HELFLAG_VISLIM_NOMOON)
     ) {
    AltM = -90; AziM = 0;
  } else {
    if (ObjectLoc(tjdut, dgeo, datm, "moon", 0, helflag, &AltM, serr) == ERR)
      return ERR;
    if (ObjectLoc(tjdut, dgeo, datm, "moon", 1, helflag, &AziM, serr) == ERR)
      return ERR;
  }
#if SWEHEL_DEBUG
{
  int i;
  for (i = 0; i < 6;i++)
    printf("dobs[%d] = %f\n", i, dobs[i]);
  printf("AltO = %.10f, AziO = %.10f\n", AltO, AziO);
  printf("AltM = %.10f, AziM = %.10f\n", AltM, AziM);
  printf("AltS = %.10f, AziS = %.10f\n", AltS, AziS);
  printf("JD = %.10f\n", tjdut);
  printf("lat = %f, eyeh = %f\n", dgeo[1], dgeo[2]);
  for (i = 0; i < 4;i++)
    printf("datm[%d] = %f\n", i, datm[i]);
  printf("helflag = %d\n", helflag);
}
#endif
  dret[0] = VisLimMagn(dobs, AltO, AziO, AltM, AziM, tjdut, AltS, AziS, sunra, dgeo[1], dgeo[2], datm, helflag, &scotopic_flag, serr);
  dret[1] = AltO;
  dret[2] = AziO;
  dret[3] = AltS;
  dret[4] = AziS;
  dret[5] = AltM;
  dret[6] = AziM;
  if (Magnitude(tjdut, dgeo, ObjectName, helflag, &(dret[7]), serr) == ERR)
    return ERR;
  retval = scotopic_flag;
  /*dret[8] = (double) is_scotopic;*/
  /*if (*serr != '\0') * in VisLimMagn(), serr is only a warning *
    retval = ERR; */
  return retval;
}

/*###################################################################
' Magn [-]
' age [Year]
' SN [-]
' AltO [deg]
' AziO [deg]
' AltM [deg]
' AziM [deg]
' MoonDistance [km]
' JDNDaysUT [-]
' AziS [deg]
' lat [deg]
' heighteye [m]
' Temperature [C]
' Pressure [mbar]
' RH [%]
' VR [km]
' TopoArcVisionis [deg]
*/
static int32 TopoArcVisionis(double Magn, double *dobs, double AltO, double AziO, double AltM, double AziM, double JDNDaysUT, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, double *dret, char *serr)
{
  double Xm, Ym, AltSi, AziSi;
  double xR = 0;
  double Xl = 45;
  double Yl, Yr;
  Yl = Magn - VisLimMagn(dobs, AltO, AziO, AltM, AziM, JDNDaysUT, AltO - Xl, AziS, sunra, Lat, HeightEye, datm, helflag, NULL, serr);
  /* if (*serr != '\0') return ERR; * serr is only a warning */
  Yr = Magn - VisLimMagn(dobs, AltO, AziO, AltM, AziM, JDNDaysUT, AltO - xR, AziS, sunra, Lat, HeightEye, datm, helflag, NULL, serr);
  /* if (*serr != '\0') return ERR; * serr is only a warning */
  /* http://en.wikipedia.org/wiki/Bisection_method*/
  if ((Yl * Yr) <= 0) {
    while(fabs(xR - Xl) > epsilon) {
      /*Calculate midpoint of domain*/
      Xm = (xR + Xl) / 2.0;
      AltSi = AltO - Xm;
      AziSi = AziS;
      Ym = Magn - VisLimMagn(dobs, AltO, AziO, AltM, AziM, JDNDaysUT, AltSi, AziSi, sunra, Lat, HeightEye, datm, helflag, NULL, serr);
      /* if (*serr != '\0') return ERR; * serr is only a warning */
      if ((Yl * Ym) > 0) {
	/* Throw away left half*/
	Xl = Xm;
	Yl = Ym;
      } else {
	/* Throw away right half */
	xR = Xm;
	Yr = Ym;
      }
    }
    Xm = (xR + Xl) / 2.0;
  } else {
    Xm = 99;
  }
  if (Xm < AltO) 
    Xm = AltO;
  *dret = Xm;
  return OK;
}

int32 CALL_CONV swe_topo_arcus_visionis(double tjdut, double *dgeo, double *datm, double *dobs, int32 helflag, double mag, double azi_obj, double alt_obj, double azi_sun, double azi_moon, double alt_moon, double *dret, char *serr)
{
  double sunra;
  swi_set_tid_acc(tjdut, helflag, 0, serr);
  sunra = SunRA(tjdut, helflag, serr);
  if (serr != NULL && *serr != '\0')
    return ERR;
  default_heliacal_parameters(datm, dgeo, dobs, helflag);
  return TopoArcVisionis(mag, dobs, alt_obj, azi_obj, alt_moon, azi_moon, tjdut, azi_sun, sunra, dgeo[1], dgeo[2], datm, helflag, dret, serr);
}

/*###################################################################*/
/*' Magn [-]
' age [Year]
' SN Snellen factor of the visual aquity of the observer
  see: http://www.i-see.org/eyecharts.html#make-your-own
' AziO [deg]
' AltM [deg]
' AziM [deg]
' MoonDistance [km]
' JDNDaysUT [-]
' AziS [deg]
' Lat [deg]
' HeightEye [m]
' Temperature [C]
' Pressure [mbar]
' RH [%]   relative humidity
' VR [km]  Meteorological Range, 
  see http://www.iol.ie/~geniet/eng/atmoastroextinction.htm
' TypeAngle 
'   [0=Object's altitude, 
'    1=Arcus Visonis (Object's altitude - Sun's altitude), 
'    2=Sun's altitude]
' HeliacalAngle [deg]
*/
static int32 HeliacalAngle(double Magn, double *dobs, double AziO, double AltM, double AziM, double JDNDaysUT, double AziS, double *dgeo, double *datm, int32 helflag, double *dangret, char *serr)
{
  double x, minx, maxx, xmin, ymin, Xl, xR, Yr, Yl, Xm, Ym, xmd, ymd;
  double Arc, DELTAx;
  double sunra = SunRA(JDNDaysUT, helflag, serr);
  double Lat = dgeo[1];
  double HeightEye = dgeo[2];
  if (PLSV == 1) {
    dangret[0] = criticalangle;
    dangret[1] = criticalangle + Magn * 2.492 + 13.447;
    dangret[2] = -(Magn * 2.492 + 13.447); /* Magn * 1.1 + 8.9;*/
    return OK;
  }
  minx = 2;
  maxx = 20;
  xmin = 0;
  ymin = 10000;
  for (x = minx; x <= maxx; x++) {
    if (TopoArcVisionis(Magn, dobs, x, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &Arc, serr) == ERR)
      return ERR;
    if (Arc < ymin) {
      ymin = Arc;
      xmin = x;
    }
  }
  Xl = xmin - 1;
  xR = xmin + 1;
  if (TopoArcVisionis(Magn, dobs, xR, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &Yr, serr) == ERR)
    return ERR;
  if (TopoArcVisionis(Magn, dobs, Xl, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &Yl, serr) == ERR)
    return ERR;
  /* http://en.wikipedia.org/wiki/Bisection_method*/
  while(fabs(xR - Xl) > 0.1) {
    /* Calculate midpoint of domain */
    Xm = (xR + Xl) / 2.0;
    DELTAx = 0.025;
    xmd = Xm + DELTAx;
    if (TopoArcVisionis(Magn, dobs, Xm, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &Ym, serr) == ERR)
      return ERR;
    if (TopoArcVisionis(Magn, dobs, xmd, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &ymd, serr) == ERR)
      return ERR;
    if (Ym >= ymd) {
      /* Throw away left half */
      Xl = Xm;
      Yl = Ym;
    } else {
      /*Throw away right half */
      xR = Xm;
      Yr = Ym;
    }
  }
  Xm = (xR + Xl) / 2.0;
  Ym = (Yr + Yl) / 2.0;
  dangret[1] = Ym;
  dangret[2] = Xm - Ym;
  dangret[0] = Xm;
  return OK;
}

int32 CALL_CONV swe_heliacal_angle(double tjdut, double *dgeo, double *datm, double *dobs, int32 helflag, double mag, double azi_obj, double azi_sun, double azi_moon, double alt_moon, double *dret, char *serr)
{
  if (dgeo[2] < SEI_ECL_GEOALT_MIN || dgeo[2] > SEI_ECL_GEOALT_MAX) {
    if (serr != NULL)
      sprintf(serr, "location for heliacal events must be between %.0f and %.0f m above sea", SEI_ECL_GEOALT_MIN, SEI_ECL_GEOALT_MAX);
    return ERR;
  }
  swi_set_tid_acc(tjdut, helflag, 0, serr);
  default_heliacal_parameters(datm, dgeo, dobs, helflag);
  return HeliacalAngle(mag, dobs, azi_obj, alt_moon, azi_moon, tjdut, azi_sun, dgeo, datm, helflag, dret, serr);
}

/*###################################################################
' AltO [deg]
' AziO [deg]
' AltS [deg]
' AziS [deg]
' parallax [deg]
' WidthMoon [deg]
*/
static double WidthMoon(double AltO, double AziO, double AltS, double AziS, double parallax)
{
  /* Yallop 1998, page 3*/
  double GeoAltO = AltO + parallax;
  return 0.27245 * parallax * (1 + sin(GeoAltO * DEGTORAD) * sin(parallax * DEGTORAD)) * (1 - cos((AltS - GeoAltO) * DEGTORAD) * cos((AziS - AziO) * DEGTORAD));
}

/*###################################################################
' W [deg]
' LengthMoon [deg]
*/
static double LengthMoon(double W, double Diamoon)
{
  double Wi, D;
  if (Diamoon == 0) Diamoon = AvgRadiusMoon * 2;
  Wi = W * 60;
  D = Diamoon * 60;
  /* Crescent length according: http://calendar.ut.ac.ir/Fa/Crescent/Data/Sultan2005.pdf*/
  return (D - 0.3 * (D + Wi) / 2.0 / Wi) / 60.0;
}

/*###################################################################
' W [deg]
' GeoARCVact [deg]
' q [-]
*/
static double qYallop(double W, double GeoARCVact)
{
  double Wi = W * 60;
  return (GeoARCVact - (11.8371 - 6.3226 * Wi + 0.7319 * Wi * Wi - 0.1018 * Wi * Wi * Wi)) / 10;
}

/*###################################################################
'A (0,p)
'B (1,q)
'C (0,r)
'D (1,s)
*/
static double crossing(double A, double B, double C, double D)
{
return (C - A) / ((B - A) - (D - C));
}

/*###################################################################*/
static int32 DeterTAV(double *dobs, double JDNDaysUT, double *dgeo, double *datm, char *ObjectName, int32 helflag, double *dret, char *serr)
{
  double Magn, AltO, AziS, AziO, AziM, AltM;
  double sunra = SunRA(JDNDaysUT, helflag, serr);
  if (Magnitude(JDNDaysUT, dgeo, ObjectName, helflag, &Magn, serr) == ERR)
    return ERR;
  if (ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 0, helflag, &AltO, serr) == ERR)
    return ERR;
  if (ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 1, helflag, &AziO, serr) == ERR)
    return ERR;
  if (strncmp(ObjectName, "moon", 4) == 0) {
    AltM = -90; 
    AziM = 0;
  } else {
    if (ObjectLoc(JDNDaysUT, dgeo, datm, "moon", 0, helflag, &AltM, serr) == ERR)
      return ERR;
    if (ObjectLoc(JDNDaysUT, dgeo, datm, "moon", 1, helflag, &AziM, serr) == ERR)
      return ERR;
  }
  if (ObjectLoc(JDNDaysUT, dgeo, datm, "sun", 1, helflag, &AziS, serr) == ERR)
    return ERR;
  if (TopoArcVisionis(Magn, dobs, AltO, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, dgeo[1], dgeo[2], datm, helflag, dret, serr) == ERR)
    return ERR;
  return OK;
}

/*###################################################################
' A y-value at x=1
' B y-value at x=0
' C y-value at x=-1
' x2min minimum for the quadratic function
*/
static double x2min(double A, double B, double C)
{
  double term = A + C - 2 * B;
  if (term == 0)
    return 0;
  return -(A - C) / 2.0 / term;
}


/*###################################################################
' A y-value at x=1
' B y-value at x=0
' C y-value at x=-1
' x
' y is y-value of quadratic function
*/
static double funct2(double A, double B, double C, double x)
{
  return (A + C - 2 * B) / 2.0 * x * x + (A - C) / 2.0 * x + B;
}

static void strcpy_VBsafe(char *sout, char *sin)
{ 
  char *sp, *sp2; 
  int iw = 0;
  sp = sin; 
  sp2 = sout;
  /* note, star name may begin with comma, such as ",zePsc" */
  while((isalnum((int) *sp) || *sp == ' ' || *sp == '-' || *sp == ',') && iw < 30) {
    *sp2 = *sp;
    sp++; sp2++; iw++;
  }
  *sp2 = '\0';
}

/*###################################################################
' JDNDaysUT [JDN]
' for information about input parameters, see function swe_heliacal_ut().
'
' output values:
'0=AltO [deg]		topocentric altitude of object (unrefracted)
'1=AppAltO [deg]        apparent altitude of object (refracted)
'2=GeoAltO [deg]        geocentric altitude of object
'3=AziO [deg]           azimuth of object
'4=AltS [deg]           topocentric altitude of Sun
'5=AziS [deg]           azimuth of Sun
'6=TAVact [deg]         actual topocentric arcus visionis
'7=ARCVact [deg]        actual (geocentric) arcus visionis
'8=DAZact [deg]         actual difference between object's and sun's azimuth
'9=ARCLact [deg]        actual longitude difference between object and sun
'10=kact [-]            extinction coefficient
'11=minTAV [deg]        smallest topocentric arcus visionis
'12=TfistVR [JDN]       first time object is visible, according to VR
'13=TbVR [JDN]          optimum time the object is visible, according to VR
'14=TlastVR [JDN]       last time object is visible, according to VR
'15=TbYallop[JDN]       best time the object is visible, according to Yallop
'16=WMoon [deg]         cresent width of moon
'17=qYal [-]            q-test value of Yallop 
'18=qCrit [-]           q-test criterion of Yallop
'19=ParO [deg]          parallax of object
'20 Magn [-]            magnitude of object
'21=RiseO [JDN]         rise/set time of object
'22=RiseS [JDN]         rise/set time of sun
'23=Lag [JDN]           rise/set time of object minus rise/set time of sun
'24=TvisVR [JDN]        visibility duration
'25=LMoon [deg]         cresent length of moon
'26=CVAact [deg]
'27=Illum [%] 'new
'28=CVAact [deg] 'new
'29=MSk [-]
*/
int32 CALL_CONV swe_heliacal_pheno_ut(double JDNDaysUT, double *dgeo, double *datm, double *dobs, char *ObjectNameIn, int32 TypeEvent, int32 helflag, double *darr, char *serr)
{
  double AziS, AltS, AltS2, AziO, AltO, AltO2, GeoAltO, AppAltO, DAZact, TAVact, ParO, MagnO;
  double ARCVact, ARCLact, kact, WMoon, LMoon = 0, qYal, qCrit;
  double RiseSetO, RiseSetS, Lag, TbYallop, TfirstVR, TlastVR, TbVR;
  double MinTAV = 0, MinTAVact, Ta, Tc, TimeStep, TimePointer, MinTAVoud = 0, DeltaAltoud = 0, DeltaAlt, TvisVR, crosspoint;
  double OldestMinTAV, extrax, illum;
  double elong, attr[30];
  double TimeCheck, LocalminCheck;
  int32 retval = OK, RS, Planet;
  AS_BOOL noriseO = FALSE;
  char ObjectName[AS_MAXCH];
  double sunra;
  int32 iflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  if (dgeo[2] < SEI_ECL_GEOALT_MIN || dgeo[2] > SEI_ECL_GEOALT_MAX) {
    if (serr != NULL)
      sprintf(serr, "location for heliacal events must be between %.0f and %.0f m above sea", SEI_ECL_GEOALT_MIN, SEI_ECL_GEOALT_MAX);
    return ERR;
  }
  swi_set_tid_acc(JDNDaysUT, helflag, 0, serr);
  sunra = SunRA(JDNDaysUT, helflag, serr);
  /* note, the fixed stars functions rewrite the star name. The input string 
     may be too short, so we have to make sure we have enough space */
  strcpy_VBsafe(ObjectName, ObjectNameIn);
  tolower_string_star(ObjectName);
  default_heliacal_parameters(datm, dgeo, dobs, helflag);
  swe_set_topo(dgeo[0], dgeo[1], dgeo[2]);
  retval = ObjectLoc(JDNDaysUT, dgeo, datm, "sun", 1, helflag, &AziS, serr);
  if (retval == OK)
    retval = ObjectLoc(JDNDaysUT, dgeo, datm, "sun", 0, helflag, &AltS, serr);
  if (retval == OK)
    retval = ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 1, helflag, &AziO, serr);
  if (retval == OK)
    retval = ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 0, helflag, &AltO, serr);
  if (retval == OK)
    retval = ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 7, helflag, &GeoAltO, serr);
  if (retval == ERR)
    return ERR;
  AppAltO = AppAltfromTopoAlt(AltO, datm[1], datm[0], helflag);
  DAZact = AziS - AziO;
  TAVact = AltO - AltS;
  /*this parallax seems to be somewhat smaller then in Yallop and SkyMap! Needs to be studied*/
  ParO = GeoAltO - AltO;
  if (Magnitude(JDNDaysUT, dgeo, ObjectName, helflag, &MagnO, serr) == ERR)
    return ERR;
  ARCVact = TAVact + ParO;
  ARCLact = acos(cos(ARCVact * DEGTORAD) * cos(DAZact * DEGTORAD)) / DEGTORAD;
  Planet = DeterObject(ObjectName);
  if (Planet == -1) {
    elong = ARCLact;
    illum = 100;
  } else {
    retval = swe_pheno_ut(JDNDaysUT, Planet, iflag|(SEFLG_TOPOCTR|SEFLG_EQUATORIAL), attr, serr);
    if (retval == ERR) return ERR;
    elong = attr[2];
    illum = attr[1] * 100;
  }
  kact = kt(AltS, sunra, dgeo[1], dgeo[2], datm[1], datm[2], datm[3], 4, serr);
  if ((0)) {
darr[26] = kR(AltS, dgeo[2]);
darr[27] = kW(dgeo[2], datm[1], datm[2]);
darr[28] = kOZ(AltS, sunra, dgeo[1]);
darr[29] = ka(AltS, sunra, dgeo[1], dgeo[2], datm[1], datm[2], datm[3], serr);
darr[30] = darr[26] + darr[27] + darr[28] + darr[29];
  }
  WMoon = 0;
  qYal = 0;
  qCrit = 0;
  LMoon = 0;
  if (Planet == SE_MOON) {
    WMoon = WidthMoon(AltO, AziO, AltS, AziS, ParO);
    LMoon = LengthMoon(WMoon, 0);
    qYal = qYallop(WMoon, ARCVact);
    if (qYal > 0.216) qCrit = 1; /* A */
    if (qYal < 0.216 && qYal > -0.014) qCrit = 2; /* B */
    if (qYal < -0.014 && qYal > -0.16) qCrit = 3; /* C */
    if (qYal < -0.16 && qYal > -0.232) qCrit = 4; /* D */
    if (qYal < -0.232 && qYal > -0.293) qCrit = 5; /* E */
    if (qYal < -0.293) qCrit = 6; /* F */
  }
  /*determine if rise or set event*/
  RS = 2;
  if (TypeEvent == 1 || TypeEvent == 4) RS = 1;
  retval = RiseSet(JDNDaysUT - 4.0 / 24.0, dgeo, datm, "sun", RS, helflag, 0, &RiseSetS, serr);
  if (retval == ERR)
    return ERR;
  retval = RiseSet(JDNDaysUT - 4.0 / 24.0, dgeo, datm, ObjectName, RS, helflag, 0, &RiseSetO, serr);
  if (retval == ERR)
    return ERR;
  TbYallop = TJD_INVALID;
  if (retval == -2) { /* object does not rise or set */
    Lag = 0;
    noriseO = TRUE;
  } else {
    Lag = RiseSetO - RiseSetS;
    if (Planet == SE_MOON)
      TbYallop = (RiseSetO * 4 + RiseSetS * 5) / 9.0;
  }
  if ((TypeEvent == 3 || TypeEvent == 4) && (Planet == -1 || Planet >= SE_MARS)) {
    TfirstVR = TJD_INVALID;
    TbVR = TJD_INVALID;
    TlastVR = TJD_INVALID;
    TvisVR = 0;
    MinTAV = 0;
    goto output_heliacal_pheno;
  }
  /* If HPheno >= 11 And HPheno <= 14 Or HPheno = 24 Then*/
  /*te bepalen m.b.v. walkthrough*/
    MinTAVact = 199;
    DeltaAlt = 0; 
    OldestMinTAV = 0;
    Ta = 0;
    Tc = 0;
    TbVR = 0;
    TimeStep = -TimeStepDefault / 24.0 / 60.0;
    if (RS == 2) TimeStep = -TimeStep;
    TimePointer = RiseSetS - TimeStep;
    do {
      TimePointer = TimePointer + TimeStep;
      OldestMinTAV = MinTAVoud;
      MinTAVoud = MinTAVact;
      DeltaAltoud = DeltaAlt;
      retval = ObjectLoc(TimePointer, dgeo, datm, "sun", 0, helflag, &AltS2, serr);
      if (retval == OK)
	retval = ObjectLoc(TimePointer, dgeo, datm, ObjectName, 0, helflag, &AltO2, serr);
      if (retval != OK)
        return ERR;
      DeltaAlt = AltO2 - AltS2;
      if (DeterTAV(dobs, TimePointer, dgeo, datm, ObjectName, helflag, &MinTAVact, serr) == ERR)
        return ERR;
      if (MinTAVoud < MinTAVact && TbVR == 0) {
        /* determine if this is a local minimum with object still above horizon*/
	TimeCheck = TimePointer + Sgn(TimeStep) * LocalMinStep / 24.0 / 60.0;
	if (RiseSetO != 0) {
	  if (TimeStep > 0)
	    TimeCheck = mymin(TimeCheck, RiseSetO);
	  else
	    TimeCheck = mymax(TimeCheck, RiseSetO);
	}
        if (DeterTAV(dobs, TimeCheck, dgeo, datm, ObjectName, helflag, &LocalminCheck, serr) == ERR)
	  return ERR;
	if (LocalminCheck > MinTAVact) {
	  extrax = x2min(MinTAVact, MinTAVoud, OldestMinTAV);
	  TbVR = TimePointer - (1 - extrax) * TimeStep;
	  MinTAV = funct2(MinTAVact, MinTAVoud, OldestMinTAV, extrax);
	}
      }
      if (DeltaAlt > MinTAVact && Tc == 0 && TbVR == 0) {
	crosspoint = crossing(DeltaAltoud, DeltaAlt, MinTAVoud, MinTAVact);
	Tc = TimePointer - TimeStep * (1 - crosspoint);
      }
      if (DeltaAlt < MinTAVact && Ta == 0 && Tc != 0) {
	crosspoint = crossing(DeltaAltoud, DeltaAlt, MinTAVoud, MinTAVact);
	Ta = TimePointer - TimeStep * (1 - crosspoint);
      }
    } while (fabs(TimePointer - RiseSetS) <= MaxTryHours / 24.0 && Ta == 0 && !((TbVR != 0 && (TypeEvent == 3 || TypeEvent == 4) && (strncmp(ObjectName, "moon", 4) != 0 && strncmp(ObjectName, "venus", 5) != 0 && strncmp(ObjectName, "mercury", 7) != 0))));
    if (RS == 2) {
      TfirstVR = Tc;
      TlastVR = Ta;
    } else {
      TfirstVR = Ta;
      TlastVR = Tc;
    }
    if (TfirstVR == 0 && TlastVR == 0) {
      if (RS == 1)
        TfirstVR = TbVR - 0.000001;
      else
        TlastVR = TbVR + 0.000001;
    }
    if (!noriseO) {
      if (RS == 1)
        TfirstVR = mymax(TfirstVR, RiseSetO);
      else
        TlastVR = mymin(TlastVR, RiseSetO);
    }
    TvisVR = TJD_INVALID; /*"#NA!" */
    if (TlastVR != 0 && TfirstVR != 0)
      TvisVR = TlastVR - TfirstVR;
    if (TlastVR == 0) TlastVR = TJD_INVALID; /*"#NA!" */
    if (TbVR == 0) TbVR = TJD_INVALID; /*"#NA!" */
    if (TfirstVR == 0) TfirstVR = TJD_INVALID; /*"#NA!" */
output_heliacal_pheno:
  /*End If*/
  darr[0] = AltO;
  darr[1] = AppAltO;
  darr[2] = GeoAltO;
  darr[3] = AziO;
  darr[4] = AltS;
  darr[5] = AziS;
  darr[6] = TAVact;
  darr[7] = ARCVact;
  darr[8] = DAZact;
  darr[9] = ARCLact;
  darr[10] = kact;
  darr[11] = MinTAV;
  darr[12] = TfirstVR;
  darr[13] = TbVR;
  darr[14] = TlastVR;
  darr[15] = TbYallop;
  darr[16] = WMoon;
  darr[17] = qYal;
  darr[18] = qCrit;
  darr[19] = ParO;
  darr[20] = MagnO;
  darr[21] = RiseSetO;
  darr[22] = RiseSetS;
  darr[23] = Lag;
  darr[24] = TvisVR;
  darr[25] = LMoon;
  darr[26] = elong;
  darr[27] = illum;
  return OK;
}

#if 0
int32 HeliacalJDut(double JDNDaysUTStart, double Age, double SN, double Lat, double Longitude, double HeightEye, double Temperature, double Pressure, double RH, double VR, char *ObjectName, int TypeEvent, char *AVkind, double *dret, char *serr)
{
  double dgeo[3], datm[4], dobs[6];
  int32 helflag = SE_HELFLAG_HIGH_PRECISION;
  helflag |= SE_HELFLAG_AVKIND_VR;
  dgeo[0] = Longitude;
  dgeo[1] = Lat;
  dgeo[2] = HeightEye;
  datm[0] = Pressure;
  datm[1] = Temperature;
  datm[2] = RH;
  datm[3] = VR;
  dobs[0] = Age;
  dobs[1] = SN;
  return swe_heliacal_ut(JDNDaysUTStart, dgeo, datm, dobs, ObjectName, TypeEvent, helflag, dret, serr);
}
#endif

static double get_synodic_period(int Planet) 
{
  /* synodic periods from:
   * Kelley/Milone/Aveni, "Exploring ancient Skies", p. 43. */
  switch(Planet) {
    case SE_MOON: return 29.530588853;
    case SE_MERCURY: return 115.8775;
    case SE_VENUS: return 583.9214;
    case SE_MARS: return 779.9361;
    case SE_JUPITER: return 398.8840;
    case SE_SATURN: return 378.0919;
    case SE_URANUS: return 369.6560;
    case SE_NEPTUNE: return 367.4867;
    case SE_PLUTO: return 366.7207;
  }
  return 366; /* for stars and default for far away planets */
}

/*###################################################################*/
static int32 moon_event_arc_vis(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, int32 TypeEvent, int32 helflag, double *dret, char *serr)
{
  double x[20], MinTAV, MinTAVoud, OldestMinTAV;
  double phase1, phase2, JDNDaysUT, JDNDaysUTi;
  double tjd_moonevent, tjd_moonevent_start;
  double DeltaAltoud, TimeCheck, LocalminCheck;
  double AltS, AltO, DeltaAlt = 90;
  char ObjectName[30];
  int32 iflag, Daystep, goingup, Planet, retval;
  int32 avkind = helflag & SE_HELFLAG_AVKIND;
  int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  dret[0] = JDNDaysUTStart; /* will be returned in error case */
  if (avkind == 0)
    avkind = SE_HELFLAG_AVKIND_VR;
  if (avkind != SE_HELFLAG_AVKIND_VR) {
    if (serr != NULL)
      strcpy(serr, "error: in valid AV kind for the moon");
    return ERR;
  }
  if (TypeEvent == 1 || TypeEvent == 2) {
    if (serr != NULL)
      strcpy(serr, "error: the moon has no morning first or evening last");
    return ERR;
  }
  strcpy(ObjectName, "moon");
  Planet = SE_MOON;
  iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
  if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
    iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
  Daystep = 1;
  if (TypeEvent == 3) {
    /*morning last */
    TypeEvent = 2;
  } else {
    /*evening first*/
    TypeEvent = 1;
    Daystep = -Daystep;
  }
  /* check Synodic/phase Period */
  JDNDaysUT = JDNDaysUTStart;
  /* start 30 days later if TypeEvent=4 (1) */
  if (TypeEvent == 1) JDNDaysUT = JDNDaysUT + 30;
  /* determination of new moon date */
  swe_pheno_ut(JDNDaysUT, Planet, iflag, x, serr);
  phase2 = x[0];
  goingup = 0;
  do {
    JDNDaysUT = JDNDaysUT + Daystep;
    phase1 = phase2;
    swe_pheno_ut(JDNDaysUT, Planet, iflag, x, serr);
    phase2 = x[0];
    if (phase2 > phase1) 
      goingup = 1;
  } while (goingup == 0 || (goingup == 1 && (phase2 > phase1)));
  /* fix the date to get the day with the smallest phase (nwest moon) */
  JDNDaysUT = JDNDaysUT - Daystep;
  /* initialize the date to look for set */
  JDNDaysUTi = JDNDaysUT;
  JDNDaysUT = JDNDaysUT - Daystep;
  MinTAVoud = 199;
  do {
    JDNDaysUT = JDNDaysUT + Daystep;
    if ((retval = RiseSet(JDNDaysUT, dgeo, datm, ObjectName, TypeEvent, helflag, 0, &tjd_moonevent, serr)) != OK)
      return retval;
    tjd_moonevent_start = tjd_moonevent;
    MinTAV = 199;
    OldestMinTAV = MinTAV;
    do {
      OldestMinTAV = MinTAVoud;
      MinTAVoud = MinTAV;
      DeltaAltoud = DeltaAlt;
      tjd_moonevent = tjd_moonevent - 1.0 / 60.0 / 24.0 * Sgn(Daystep);
      if (ObjectLoc(tjd_moonevent, dgeo, datm, "sun", 0, helflag, &AltS, serr) == ERR)
	return ERR;
      if (ObjectLoc(tjd_moonevent, dgeo, datm, ObjectName, 0, helflag, &AltO, serr) == ERR)
	return ERR;
      DeltaAlt = AltO - AltS;
      if (DeterTAV(dobs, tjd_moonevent, dgeo, datm, ObjectName, helflag, &MinTAV, serr) == ERR)
        return ERR;
      TimeCheck = tjd_moonevent - LocalMinStep / 60.0 / 24.0 * Sgn(Daystep);
      if (DeterTAV(dobs, TimeCheck, dgeo, datm, ObjectName, helflag, &LocalminCheck, serr) == ERR)
        return ERR;
/*printf("%f, %f <= %f\n", tjd_moonevent, MinTAV, MinTAVoud);*/
    /* while (MinTAV <= MinTAVoud && fabs(tjd_moonevent - tjd_moonevent_start) < 120.0 / 60.0 / 24.0);*/
    } while ((MinTAV <= MinTAVoud || LocalminCheck < MinTAV) && fabs(tjd_moonevent - tjd_moonevent_start) < 120.0 / 60.0 / 24.0);
  /* while (DeltaAlt < MinTAVoud && fabs(JDNDaysUT - JDNDaysUTi) < 15);*/
  } while (DeltaAltoud < MinTAVoud && fabs(JDNDaysUT - JDNDaysUTi) < 15);
  if (fabs(JDNDaysUT - JDNDaysUTi) < 15) {
    tjd_moonevent += (1 - x2min(MinTAV, MinTAVoud, OldestMinTAV)) * Sgn(Daystep) / 60.0 / 24.0;
  } else {
    strcpy(serr, "no date found for lunar event");
    return ERR;
  }
  dret[0] = tjd_moonevent;
  return OK;
}

static int32 heliacal_ut_arc_vis(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 TypeEventIn, int32 helflag, double *dret, char *serr_ret)
{
  double x[6];
  double xin[2];
  double xaz[2];
  double dang[3];
  double objectmagn = 0, maxlength, DayStep;
  double JDNDaysUT, JDNDaysUTfinal, JDNDaysUTstep, JDNDaysUTstepoud, JDNarcvisUT, tjd_tt, tret, OudeDatum, JDNDaysUTinp = JDNDaysUTStart, JDNDaysUTtijd;
  double ArcusVis, ArcusVisDelta, ArcusVisPto, ArcusVisDeltaoud;
  double Trise, sunsangle, Theliacal, Tdelta, Angle;
  double TimeStep, TimePointer, OldestMinTAV, MinTAVoud, MinTAVact, extrax, TbVR = 0;
  double AziS, AltS, AziO, AltO, DeltaAlt;
  double direct, Pressure, Temperature, d;
  int32 epheflag, retval = OK;
  int32 iflag, Planet, eventtype;
  int32 TypeEvent = TypeEventIn;
  int doneoneday;
  char serr[AS_MAXCH];
  *dret = JDNDaysUTStart;
  *serr = '\0';
  Planet = DeterObject(ObjectName);
  Pressure = datm[0];
  Temperature = datm[1];
  /* determine Magnitude of star*/
  if ((retval = Magnitude(JDNDaysUTStart, dgeo, ObjectName, helflag, &objectmagn, serr)) == ERR)
    goto swe_heliacal_err;
  epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
  if (!(helflag & SE_HELFLAG_HIGH_PRECISION)) 
    iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
  /* start values for search of heliacal rise
   * maxlength = phase period in days, smaller than minimal synodic period */
  /* days per step (for heliacal rise) in power of two */
  switch(Planet) {
    case SE_MERCURY: 
      DayStep = 1; maxlength = 100; break;
    case SE_VENUS: 
      DayStep = 64; maxlength = 384; break;
    case SE_MARS: 
      DayStep = 128; maxlength = 640; break;
    case SE_JUPITER: 
      DayStep = 64; maxlength = 384; break;
    case SE_SATURN: 
      DayStep = 64; maxlength = 256; break;
    default:
      DayStep = 64; maxlength = 256; break;
  }
  /* heliacal setting */
  eventtype = TypeEvent;
  if (eventtype == 2) DayStep = -DayStep;
  /* acronychal setting */
  if (eventtype == 4) {
      eventtype = 1;
      DayStep = -DayStep;
  }
  /* acronychal rising */
  if (eventtype == 3) eventtype = 2;
  eventtype |= SE_BIT_DISC_CENTER;
  /* normalize the maxlength to the step size */
  {
    /* check each Synodic/phase Period */
    JDNDaysUT = JDNDaysUTStart;
    /* make sure one can find an event on the just after the JDNDaysUTStart */
    JDNDaysUTfinal = JDNDaysUT + maxlength;
    JDNDaysUT = JDNDaysUT - 1;
    if (DayStep < 0) {
      JDNDaysUTtijd = JDNDaysUT;
      JDNDaysUT = JDNDaysUTfinal;
      JDNDaysUTfinal = JDNDaysUTtijd;
    }
    /* prepair the search */
    JDNDaysUTstep = JDNDaysUT - DayStep;
    doneoneday = 0;
    ArcusVisDelta = 199;
    ArcusVisPto = -5.55;
    do { /* this is a do {} while() loop */
      if (fabs(DayStep) == 1) doneoneday = 1;
      do { /* this is a do {} while() loop */
	/* init search for heliacal rise */
	JDNDaysUTstepoud = JDNDaysUTstep;
	ArcusVisDeltaoud = ArcusVisDelta;
	JDNDaysUTstep = JDNDaysUTstep + DayStep;
	/* determine rise/set time */
	if ((retval = my_rise_trans(JDNDaysUTstep, SE_SUN, "", eventtype, helflag, dgeo, datm, &tret, serr)) == ERR)
	  goto swe_heliacal_err;
	/* determine time compensation to get Sun's altitude at heliacal rise */
	tjd_tt = tret + swe_deltat_ex(tret, epheflag, serr);
	if ((retval = swe_calc(tjd_tt, SE_SUN, iflag, x, serr)) == ERR)
	  goto swe_heliacal_err;
	xin[0] = x[0];
	xin[1] = x[1];
	swe_azalt(tret, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
	Trise = HourAngle(xaz[1], x[1], dgeo[1]);
	sunsangle = ArcusVisPto;
	if (helflag & SE_HELFLAG_AVKIND_MIN7) sunsangle = -7;
	if (helflag & SE_HELFLAG_AVKIND_MIN9) sunsangle = -9;
	Theliacal = HourAngle(sunsangle, x[1], dgeo[1]);
	Tdelta = Theliacal - Trise;
	if (TypeEvent == 2 || TypeEvent== 3) Tdelta = -Tdelta;
	/* determine appr.time when sun is at the wanted Sun's altitude */
	JDNarcvisUT = tret - Tdelta / 24;
	tjd_tt = JDNarcvisUT + swe_deltat_ex(JDNarcvisUT, epheflag, serr);
	/* determine Sun's position */
	if ((retval = swe_calc(tjd_tt, SE_SUN, iflag, x, serr)) == ERR)
	  goto swe_heliacal_err;
	xin[0] = x[0];
	xin[1] = x[1];
	swe_azalt(JDNarcvisUT, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
	AziS = xaz[0] + 180;
	if (AziS >= 360) AziS = AziS - 360;
	AltS = xaz[1];
	/* determine Moon's position */
#if 0
  double AltM, AziM;
	if ((retval = swe_calc(tjd_tt, SE_MOON, iflag, x, serr)) == ERR)
	  goto swe_heliacal_err;
	xin[0] = x[0];
	xin[1] = x[1];
	swe_azalt(JDNarcvisUT, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
	AziM = xaz[0] + 180;
	if (AziM >= 360) AziM = AziM - 360;
	AltM = xaz[1];
#endif
	/* determine object's position */
	if (Planet != -1) {
	  if ((retval = swe_calc(tjd_tt, Planet, iflag, x, serr)) == ERR)
	    goto swe_heliacal_err;
	  /* determine magnitude of Planet */
	  if ((retval = Magnitude(JDNarcvisUT, dgeo, ObjectName, helflag, &objectmagn, serr)) == ERR)
	    goto swe_heliacal_err;
	} else {
	  if ((retval = call_swe_fixstar(ObjectName, tjd_tt, iflag, x, serr)) == ERR)
	    goto swe_heliacal_err;
	}
	xin[0] = x[0];
	xin[1] = x[1];
	swe_azalt(JDNarcvisUT, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
	AziO = xaz[0] + 180;
	if (AziO >= 360) AziO = AziO - 360;
	AltO = xaz[1];
	/* determine arcusvisionis */
	DeltaAlt = AltO - AltS;
	/*if ((retval = HeliacalAngle(objectmagn, dobs, AziO, AltM, AziM, JDNarcvisUT, AziS, dgeo, datm, helflag, dang, serr)) == ERR)*/
	if ((retval = HeliacalAngle(objectmagn, dobs, AziO, -1, 0, JDNarcvisUT, AziS, dgeo, datm, helflag, dang, serr)) == ERR)
	  goto swe_heliacal_err;
	ArcusVis = dang[1];
	ArcusVisPto = dang[2];
	ArcusVisDelta = DeltaAlt - ArcusVis;
      /*} while (((ArcusVisDeltaoud > 0 && ArcusVisDelta < 0) || ArcusVisDelta < 0) && (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep) > 0);*/
      } while ((ArcusVisDeltaoud > 0 || ArcusVisDelta < 0) && (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep) > 0);
      if (doneoneday == 0 && (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep) > 0) {
	/* go back to date before heliacal altitude */
	ArcusVisDelta = ArcusVisDeltaoud;
	DayStep = ((int) (fabs(DayStep) / 2.0)) * Sgn(DayStep);
	JDNDaysUTstep = JDNDaysUTstepoud;
      }
    } while (doneoneday == 0 && (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep) > 0);
  }
  d = (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep);
  if (d <= 0 || d >= maxlength) {
    dret[0] = JDNDaysUTinp; /* no date found, just return input */
    retval = -2; /* marks "not found" within synodic period */
    sprintf(serr, "heliacal event not found within maxlength %f\n", maxlength);
    goto swe_heliacal_err;
  } 
#if 0
  if (helflag & SE_HELFLAG_AVKIND_VR) {
    double darr[40];
    if (swe_heliacal_pheno_ut(JDNarcvisUT, dgeo, datm, dobs, ObjectName, TypeEvent, helflag, darr, serr) != OK)
      return ERR;
    JDNarcvisUT = darr[13];
    }
  }
#endif
  direct = TimeStepDefault / 24.0 / 60.0;
  if (DayStep < 0) direct = -direct;
  if (helflag & SE_HELFLAG_AVKIND_VR) {
    /*te bepalen m.b.v. walkthrough*/
    TimeStep = direct;
    TbVR = 0;
    TimePointer = JDNarcvisUT;
    if (DeterTAV(dobs, TimePointer, dgeo, datm, ObjectName, helflag, &OldestMinTAV, serr) == ERR)
      return ERR;
    TimePointer = TimePointer + TimeStep;
    if (DeterTAV(dobs, TimePointer, dgeo, datm, ObjectName, helflag, &MinTAVoud, serr) == ERR)
      return ERR;
    if (MinTAVoud > OldestMinTAV) {
      TimePointer = JDNarcvisUT;
      TimeStep = -TimeStep;
      MinTAVact = OldestMinTAV;
    } else {
      MinTAVact = MinTAVoud;
      MinTAVoud = OldestMinTAV;
    }
    /*TimePointer = TimePointer - Timestep*/
    do {
      TimePointer = TimePointer + TimeStep;
      OldestMinTAV = MinTAVoud;
      MinTAVoud = MinTAVact;
      if (DeterTAV(dobs, TimePointer, dgeo, datm, ObjectName, helflag, &MinTAVact, serr) == ERR)
        return ERR;
      if (MinTAVoud < MinTAVact) {
	extrax = x2min(MinTAVact, MinTAVoud, OldestMinTAV);
	TbVR = TimePointer - (1 - extrax) * TimeStep;
      }
    } while (TbVR == 0);
    JDNarcvisUT = TbVR;
  }
  /*if (strncmp(AVkind, "pto", 3) == 0) */
  if (helflag & SE_HELFLAG_AVKIND_PTO) {
    do {
      OudeDatum = JDNarcvisUT;
      JDNarcvisUT = JDNarcvisUT - direct;
      tjd_tt = JDNarcvisUT + swe_deltat_ex(JDNarcvisUT, epheflag, serr);
      if (Planet != -1) {
	if ((retval = swe_calc(tjd_tt, Planet, iflag, x, serr)) == ERR)
	  goto swe_heliacal_err;
      } else {
	if ((retval = call_swe_fixstar(ObjectName, tjd_tt, iflag, x, serr)) == ERR)
	  goto swe_heliacal_err;
      }
      xin[0] = x[0];
      xin[1] = x[1];
      swe_azalt(JDNarcvisUT, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
      Angle = xaz[1];
    } while (Angle > 0);
    JDNarcvisUT = (JDNarcvisUT + OudeDatum) / 2.0;
  }
  if (JDNarcvisUT < -9999999 || JDNarcvisUT > 9999999) {
    dret[0] = JDNDaysUT; /* no date found, just return input */
    strcpy(serr, "no heliacal date found");
    retval = ERR;
    goto swe_heliacal_err;
  }
  dret[0] = JDNarcvisUT;
swe_heliacal_err:
  if (serr_ret != NULL && *serr != '\0')
    strcpy(serr_ret, serr);
  return retval;
}

static int32 get_asc_obl(double tjd, int32 ipl, char *star, int32 iflag, double
 *dgeo, AS_BOOL desc_obl, double *daop, char *serr)
{
  int32 retval;
  int32 epheflag = iflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  double x[6], adp;
  char s[AS_MAXCH];
  char star2[AS_MAXCH];
  strcpy(star2, star);
  if (ipl == -1) {
    if ((retval = swe_fixstar(star2, tjd, epheflag | SEFLG_EQUATORIAL, x, serr)) == ERR)
      return ERR;
  } else {
    if ((retval = swe_calc(tjd, ipl, epheflag | SEFLG_EQUATORIAL, x, serr)) == ERR)
      return ERR;
  }
  adp = tan(dgeo[1] * DEGTORAD) * tan(x[1] * DEGTORAD);
  if (fabs(adp) > 1) {
    if (star != NULL && *star != '\0')
      strcpy(s, star);
    else
      swe_get_planet_name(ipl, s);
    sprintf(serr, "%s is circumpolar, cannot calculate heliacal event", s);
    return -2;
  }
  adp = asin(adp) / DEGTORAD;
  if (desc_obl)
    *daop = x[0] + adp;
  else
    *daop = x[0] - adp;
  *daop = swe_degnorm(*daop);
  return OK;
}

#if 0
static int32 get_asc_obl_old(double tjd, int32 ipl, char *star, int32 iflag, double *dgeo, AS_BOOL desc_obl, double *daop, char *serr)
{
  int32 retval;
  int32 epheflag = iflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  double x[6], adp;
  char s[AS_MAXCH];
  if (star != NULL && *star != '\0') {
    if ((retval = call_swe_fixstar(star, tjd, epheflag | SEFLG_EQUATORIAL, x, serr)) == ERR)
      return ERR;
  } else {
    if ((retval = swe_calc(tjd, ipl, epheflag | SEFLG_EQUATORIAL, x, serr)) == ERR)
      return ERR;
  }
  adp = tan(dgeo[1] * DEGTORAD) * tan(x[1] * DEGTORAD); 
  if (fabs(adp) > 1) {
    if (star != NULL && *star != '\0') 
      strcpy(s, star);
    else 
      swe_get_planet_name(ipl, s);
    sprintf(serr, "%s is circumpolar, cannot calculate heliacal event", s);
    return -2;
  }
  adp = asin(adp) / DEGTORAD;
  if (desc_obl)
    *daop = x[0] + adp;
  else
    *daop = x[0] - adp;
  *daop = swe_degnorm(*daop);
  return OK;
}
#endif

static int32 get_asc_obl_diff(double tjd, int32 ipl, char *star, int32 iflag, double *dgeo, AS_BOOL desc_obl, AS_BOOL is_acronychal, double *dsunpl, char *serr)
{
  int32 retval = OK;
  double aosun, aopl;
  /* ascensio obliqua of sun */
  retval = get_asc_obl(tjd, SE_SUN, "", iflag, dgeo, desc_obl, &aosun, serr);
  if (retval != OK)
    return retval;
  if (is_acronychal) {
    if (desc_obl == TRUE)
      desc_obl = FALSE;
    else
      desc_obl = TRUE;
  }
  /* ascensio obliqua of body */
  retval = get_asc_obl(tjd, ipl, star, iflag, dgeo, desc_obl, &aopl, serr);
  if (retval != OK)
    return retval;
  *dsunpl = swe_degnorm(aosun - aopl);
  if (is_acronychal)
    *dsunpl = swe_degnorm(*dsunpl - 180);
  if (*dsunpl > 180) *dsunpl -= 360;
  return OK;
}

#if 0
static int32 get_asc_obl_diff_old(double tjd, int32 ipl, char *star, int32 iflag, double *dgeo, AS_BOOL desc_obl, double *dsunpl, char *serr)
{
  int32 retval = OK;
  double aosun, aopl;
  /* ascensio obliqua of sun */
  retval = get_asc_obl(tjd, SE_SUN, "", iflag, dgeo, desc_obl, &aosun, serr);
  if (retval != OK)
    return retval;
  /* ascensio obliqua of body */
  retval = get_asc_obl(tjd, ipl, star, iflag, dgeo, desc_obl, &aopl, serr);
  if (retval != OK)
    return retval;
  *dsunpl = swe_degnorm(aosun - aopl);
  return OK;
}
#endif

/* times of 
 * - superior and inferior conjunction (Mercury and Venus)
 * - conjunction and opposition (ipl >= Mars)
 */
static const double tcon[] =
{
  0, 0, 
  2451550, 2451550,  /* Moon */
  2451604, 2451670,  /* Mercury */
  2451980, 2452280,  /* Venus */
  2451727, 2452074,  /* Mars */
  2451673, 2451877,  /* Jupiter */ 
  2451675, 2451868,  /* Saturn */
  2451581, 2451768,  /* Uranus */ 
  2451568, 2451753,  /* Neptune */
};

static int32 find_conjunct_sun(double tjd_start, int32 ipl, int32 helflag, int32 TypeEvent, double *tjd, char *serr)
{
  int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  int i;
  double tjdcon, tjd0, ds, dsynperiod, x[6], xs[6], daspect = 0;
  if (ipl >= SE_MARS && TypeEvent >= 3)
    daspect = 180;
  i = (TypeEvent - 1) / 2 + ipl * 2;
  tjd0 = tcon[i];
  dsynperiod = get_synodic_period(ipl);
  tjdcon = tjd0 + ((floor) ((tjd_start - tjd0) / dsynperiod) + 1) * dsynperiod;
  ds = 100;
  while (ds > 0.5) {
    if (swe_calc(tjdcon, ipl, epheflag|SEFLG_SPEED, x, serr) == ERR)
      return ERR;
    if (swe_calc(tjdcon, SE_SUN, epheflag|SEFLG_SPEED, xs, serr) == ERR)
      return ERR;
    ds = swe_degnorm(x[0] - xs[0] - daspect);
    if (ds > 180) ds -= 360;
    tjdcon -= ds / (x[3] - xs[3]);
  }
  *tjd = tjdcon;
  return OK;
}

static int32 get_asc_obl_with_sun(double tjd_start, int32 ipl, char *star, int32 helflag, int32 evtyp, double dperiod, double *dgeo, double *tjdret, char *serr)
{
  int i, retval;
  int32 is_acronychal = FALSE;
  int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  double dsunpl = 1, dsunpl_save, dsunpl_test, tjd, daystep;
  AS_BOOL desc_obl = FALSE, retro = FALSE;
  if (evtyp == SE_EVENING_LAST || evtyp == SE_EVENING_FIRST)
    desc_obl = TRUE;
  if (evtyp == SE_MORNING_FIRST || evtyp == SE_EVENING_LAST)
    retro = TRUE;
  if (evtyp == SE_ACRONYCHAL_RISING)
    desc_obl = TRUE;
  if (evtyp == SE_ACRONYCHAL_RISING || evtyp ==  SE_ACRONYCHAL_SETTING) {
    is_acronychal = TRUE;
    if (ipl != SE_MOON)
      retro = TRUE;
  }
  //  if (evtyp == 3 || evtyp == 4)
  //    dangsearch = 180;
  /* find date when sun and object have the same ascensio obliqua */
  tjd = tjd_start;
  dsunpl_save = -999999999;
  retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, is_acronychal, &dsunpl, serr);
  if (retval != OK)  /* retval may be ERR or -2 */
    return retval;
  daystep = 20;
  i = 0;
  while (dsunpl_save == -999999999 ||
      /*fabs(dsunpl - dsunpl_save) > 180 ||*/
      fabs(dsunpl) + fabs(dsunpl_save) > 180 ||
      (retro && !(dsunpl_save < 0 && dsunpl >= 0)) ||
      (!retro && !(dsunpl_save >= 0 && dsunpl < 0))) {
    i++;
    if (i > 5000) {
      sprintf(serr, "loop in get_asc_obl_with_sun() (1)");
      return ERR;
    }
    dsunpl_save = dsunpl;
    tjd += 10.0;
    if (dperiod > 0 && tjd - tjd_start > dperiod)
      return -2;
    retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, is_acronychal, &dsunpl, serr);
    if (retval != OK)  /* retval may be ERR or -2 */
      return retval;
  }
  tjd_start = tjd - daystep;
  daystep /= 2.0;
  tjd = tjd_start + daystep;
  retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, is_acronychal, &dsunpl_test, serr);
  if (retval != OK)  /* retval may be ERR or -2 */
    return retval;
  i = 0;
  while (fabs(dsunpl) > 0.00001) {
    i++;
    if (i > 5000) {
      sprintf(serr, "loop in get_asc_obl_with_sun() (2)");
      return ERR;
    }
    if (dsunpl_save * dsunpl_test >= 0) {
      dsunpl_save = dsunpl_test;
      tjd_start = tjd;
    } else {
      dsunpl = dsunpl_test;
    }
    daystep /= 2.0;
    tjd = tjd_start + daystep;
    retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, is_acronychal, &dsunpl_test, serr);
    if (retval != OK)  /* retval may be ERR or -2 */
      return retval;
  }
  *tjdret = tjd;
  return OK;
}

#if 0
/* works only for fixed stars */
static int32 get_asc_obl_with_sun_old(double tjd_start, int32 ipl, char *star, int32 helflag, int32 TypeEvent, double *dgeo, double *tjdret, char *serr)
{
  int retval;
  int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  double dsunpl = 1, tjd, daystep, dsunpl_save;
  double dsynperiod = 367;
  double dangsearch = 0;
  AS_BOOL desc_obl = FALSE;
  if (TypeEvent == 2 || TypeEvent == 3)
    desc_obl = TRUE;
  if (TypeEvent == 3 || TypeEvent == 4)
    dangsearch = 180;
  /* find date when sun and object have the same ascensio obliqua */
  daystep = dsynperiod;
  tjd = tjd_start;
  retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, &dsunpl, serr);
  if (retval != OK)  /* retval may be ERR or -2 */
    return retval;
  while (dsunpl < 359.99999) {
    dsunpl_save = dsunpl;
    daystep /= 2.0;
    retval = get_asc_obl_diff(tjd + daystep, ipl, star, epheflag, dgeo, desc_obl, &dsunpl, serr);
    if (retval != OK)  /* retval may be ERR or -2 */
      return retval;
    if (dsunpl > dsunpl_save)
      tjd += daystep;
    else
      dsunpl = dsunpl_save;
  }
  *tjdret = tjd;
  return OK;
}
#endif

#if 0
/* works only for fixed stars */
static int32 get_asc_obl_acronychal(double tjd_start, int32 ipl, char *star, int32 helflag, int32 TypeEvent, double *dgeo, double *tjdret, char *serr)
{
  int retval;
  int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  double dsunpl = 1, tjd, daystep, dsunpl_save;
  double dsynperiod = 367;
  double aosun, aopl;
  AS_BOOL sun_desc = TRUE, obj_desc = FALSE;
  daystep = dsynperiod;
  tjd = tjd_start;
  if (TypeEvent == 4) {
    sun_desc = FALSE; 
    obj_desc = TRUE;
  }
  /* ascensio (descensio) obliqua of sun */
  retval = get_asc_obl(tjd, SE_SUN, "", epheflag, dgeo, sun_desc, &aosun, serr);
  if (retval != OK)  /* retval may be ERR or -2 */
    return retval;
  /* ascensio (descensio) obliqua of body */
  retval = get_asc_obl(tjd, ipl, star, epheflag, dgeo, obj_desc, &aopl, serr);
  if (retval != OK)  /* retval may be ERR or -2 */
    return retval;
  dsunpl = swe_degnorm(aosun - aopl + 180);
  while (dsunpl < 359.99999) {
    dsunpl_save = dsunpl;
    daystep /= 2.0;
    /* ascensio (descensio) obliqua of sun */
    retval = get_asc_obl(tjd+daystep, SE_SUN, "", epheflag, dgeo, sun_desc, &aosun, serr);
    if (retval != OK)  /* retval may be ERR or -2 */
      return retval;
    /* ascensio (descensio) obliqua of body */
    retval = get_asc_obl(tjd+daystep, ipl, star, epheflag, dgeo, obj_desc, &aopl, serr);
    if (retval != OK)  /* retval may be ERR or -2 */
      return retval;
    dsunpl = swe_degnorm(aosun - aopl + 180);
    if (dsunpl > dsunpl_save)
      tjd += daystep;
    else
      dsunpl = dsunpl_save;
  }
  *tjdret = tjd;
  return OK;
}
#endif

static int32 get_heliacal_day(double tjd, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, int32 TypeEvent, double *thel, char *serr)
{
  int32 i, visible_at_sunsetrise, is_rise_or_set = 0, ndays, retval, retval_old;
  double direct_day = 0, direct_time = 0, tfac, tend, daystep, tday, vdelta, tret;
  double darr[30], vd, dmag, div;
  int32 ipl = DeterObject(ObjectName);
  /* 
   * find the day and minute on which the object becomes visible 
   */
  switch (TypeEvent) {
    /* morning first */
    case 1: is_rise_or_set = SE_CALC_RISE; 
      direct_day = 1; direct_time = -1;
      break;
    /* evening last */
    case 2: is_rise_or_set = SE_CALC_SET; 
      direct_day = -1; direct_time = 1;
      break;
    /* evening first */
    case 3: is_rise_or_set = SE_CALC_SET; 
      direct_day = 1; direct_time = 1;
      break;
    /* morning last */
    case 4: is_rise_or_set = SE_CALC_RISE; 
      direct_day = -1; direct_time = -1;
      break;
  }
  tfac = 1;
  switch (ipl) {
    case SE_MOON: 
      ndays = 16; 
      daystep = 1;
      break;
    case SE_MERCURY: 
      ndays = 60; tjd -= 0 * direct_day; 
      daystep = 5;
      tfac = 5;
      break;
    case SE_VENUS: 
      ndays = 300; tjd -= 30 * direct_day; 
      daystep = 5;
      if (TypeEvent >= 3) {
	daystep = 15;
        tfac = 3;
      }
      break;
    case SE_MARS: 
      ndays = 400; 
      daystep = 15;
      tfac = 5;
      break; 
    case SE_SATURN: 
      ndays = 300; 
      daystep = 20;
      tfac = 5;
      break; 
    case -1:
      ndays = 300;
      if (call_swe_fixstar_mag(ObjectName, &dmag, serr) == ERR) {
	return ERR;
      }
      daystep = 15;
      tfac = 10;
      if (dmag > 2) {
        daystep = 15;
      }
      if (dmag < 0) {
	tfac = 3;
      }
      break;
    default:
      ndays = 300; 
      daystep = 15;
      tfac = 3;
      break;
  }
  tend = tjd + ndays * direct_day;
  retval_old = -2;
  for (tday = tjd, i = 0;
       (direct_day > 0 && tday < tend) || (direct_day < 0 && tday > tend);
       tday += daystep * direct_day, i++) {
    vdelta = -100; 
    if (i > 0)
      tday -= 0.3 * direct_day;
    if ((retval = my_rise_trans(tday, SE_SUN, "", is_rise_or_set, helflag, dgeo, datm, &tret, serr)) == ERR) {
      return ERR;
    }
    /* sun does not rise: try next day */
    if (retval == -2) {
      retval_old = retval;
      continue;
    }
    retval = swe_vis_limit_mag(tret, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
    if (retval == ERR)
      return ERR;
#if 1
    /*  object has appeared above horizon: reduce daystep */
    if (retval_old == -2 && retval >= 0 && daystep > 1) {
      retval_old = retval;
      tday -= daystep * direct_day;
      daystep = 1;
      /* Note: beyond latitude 55N (?), Mars can have a morning last. 
       * If the period of visibility is less than 5 days, we may miss the
       * event. I don't know if this happens */
      if (ipl >= SE_MARS || ipl == -1)
        daystep = 5;
      continue;
    }
    retval_old = retval;
#endif
    /*  object below horizon: try next day */
    if (retval == -2)
      continue;
    vdelta = darr[0] - darr[7];
    /* find minute of object's becoming visible */
    div = 1440.0;
//    div = 86400.0;
//      printf("tret=%f, m1=%f, m2=%f\n", tret, darr[0], darr[1]);
    vd = -1;
    visible_at_sunsetrise = 1;
    while (retval != -2 && (vd = darr[0] - darr[7]) < 0) {
      visible_at_sunsetrise = 0;
      if (vd < -1.0)
	tret += 5.0 / div * direct_time * tfac;
      else if (vd < -0.5)
	tret += 2.0 / div * direct_time * tfac;
      else if (vd < -0.1)
	tret += 1.0 / div * direct_time * tfac;
      else 
	tret += 1.0 / div * direct_time;
      retval = swe_vis_limit_mag(tret, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
      if (retval == ERR)
	return ERR;
    }
    /* if possible move a bit away from sunset, where vis_limit_mag() has strange behaviour */
    if (visible_at_sunsetrise) {
      for (i = 0; i < 10; i++) {
	if ((retval = swe_vis_limit_mag(tret + 1.0 / div * direct_time, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) >= 0
	  && darr[0] - darr[7] > vd) {
          vd = darr[0] - darr[7];
	  tret += 1.0 / div * direct_time;
	}
      }
    }
    vdelta = darr[0] - darr[7];
    /* object is visible, save time of appearance */
    if (vdelta > 0) {
      if ((ipl >= SE_MARS || ipl == -1) && daystep > 1) {
	tday -= daystep * direct_day;
	daystep = 1;
      } else {
	*thel = tret;
	return OK;
      }
    }
  }
  sprintf(serr, "heliacal event does not happen");
  //printf("%s\n", serr);
  return -2;
}

static int32 time_optimum_visibility(double tjd, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, double *tret, char *serr)
{
  int32 retval, retval_sv, i;
  double t1, t2, vl1, vl2, d, darr[10], phot_scot_opic, phot_scot_opic_sv;
  // double vl;
  int t_has_changed;
  *tret = tjd;
  retval = swe_vis_limit_mag(tjd, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
  if (retval == ERR) return ERR;
  retval_sv = retval;
  //vl = darr[0] - darr[7];
  //vl = -1;
  t1 = tjd;
  t2 = tjd;
  vl1 = -1; 
  vl2 = -1;
  //printf("begin tret=%f, dvl=%f\n", tjd, darr[0] - darr[7]);
  phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
  for (i = 0, d = 100.0 / 86400.0; i < 3; i++, d /= 10.0) {
    // fprintf(stderr, "i= %d\n", i);
    t1 += d;
    t_has_changed = 0;
    while((retval = swe_vis_limit_mag(t1 - d, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) >= 0 
	&& darr[0] > darr[7] 
	&& darr[0] - darr[7] > vl1) { 
      t1 -= d; vl1 = darr[0] - darr[7]; 
      t_has_changed = 1;
      //fprintf(stderr, "vl1=%f %d vlm=%f, obm=%f, t=%f\n", vl, retval, darr[0], darr[7], tjd + d);
      retval_sv = retval;
      phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
    /*  printf("1: %f\n", darr[8]);*/
    }
    if (t_has_changed == 0)
      t1 -= d; /* revert initial addition */
    if (retval == ERR) return ERR;
  }
  for (i = 0, d = 100.0 / 86400.0; i < 3; i++, d /= 10.0) {
    t2 -= d;
    t_has_changed = 0;
    while((retval = swe_vis_limit_mag(t2 + d, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) >= 0 
        && darr[0] > darr[7] 
	&& darr[0] - darr[7] > vl2) { 
      t2 += d; vl2 = darr[0] - darr[7]; 
      t_has_changed = 1;
      //fprintf(stderr, "vl2=%f %d vlm=%f, obm=%f, t=%f\n", vl, retval, darr[0], darr[7], tjd + d);
      retval_sv = retval;
      phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
    /*  printf("2: %f\n", darr[8]);*/
    }
    if (t_has_changed == 0)
      t2 += d; /* revert initial subtraction */
    if (retval == ERR) return ERR;
  }
  if (vl2 > vl1) 
    tjd = t2;
  else
    tjd = t1;
    /*  printf("3: %f <-> %f\n", darr[8], phot_scot_opic_sv);*/
  *tret = tjd;
  //printf("end tret=%f, dvl=%f\n", tjd, darr[0] - darr[7]);
  if (retval >= 0) {
    /* search for optimum came to an end because change scotopic/photopic: */
    phot_scot_opic = (retval & SE_SCOTOPIC_FLAG);
    if (phot_scot_opic_sv != phot_scot_opic) {
      /* calling function writes warning into serr */
      printf ("hallo -2\n");
      return -2;
    }
    /* valid result found but it is close to the scotopic/photopic limit */
    if (retval_sv & SE_MIXEDOPIC_FLAG) {
      printf ("hallo -2\n");
      return -2;
    }
  }
  return OK;
}

static int32 time_limit_invisible(double tjd, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, int32 direct, double *tret, char *serr)
{
  int32 retval, retval_sv, i, ncnt = 3;
  double d = 0, darr[10], phot_scot_opic, phot_scot_opic_sv;
  double d0 = 100.0 / 86400.0;
  *tret = tjd;
  if (strcmp(ObjectName, "moon") == 0) {
    d0 *= 10;
    ncnt = 4;
  }
  retval = swe_vis_limit_mag(tjd + d * direct, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
  if (retval == ERR) return ERR;
  retval_sv = retval;
  phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
  for (i = 0, d = d0; i < ncnt; i++, d /= 10.0) {
    while((retval = swe_vis_limit_mag(tjd + d * direct, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) >= 0 
        && darr[0] > darr[7]) { 
      tjd += d * direct; 
      retval_sv = retval;
      phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
  /*    printf("%d: %f\n", direct, darr[8]); */
    }
  }
   /*   printf("4: %f, %f/%f %f <-> %f\n", darr[8], darr[0], darr[7], tjd, phot_scot_opic_sv); */
  *tret = tjd;
  /* if object disappears at setting, retval is -2, but we want it OK, and
   * also suppress the warning "object is below local horizon" */
  *serr = '\0'; 
  if (retval >= 0) {
    /* search for limit came to an end because change scotopic/photopic: */
    phot_scot_opic = (retval & SE_SCOTOPIC_FLAG);
    if (phot_scot_opic_sv != phot_scot_opic) {
      /* calling function writes warning into serr */
      return -2;
    }
    /* valid result found but it is close to the scotopic/photopic limit */
    if (retval_sv & SE_MIXEDOPIC_FLAG) {
      return -2;
    }
  }
  return OK;
}

static int32 get_acronychal_day(double tjd, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, int32 TypeEvent, double *thel, char *serr) {
  double tret, tret_dark, darr[30], dtret;
  /* x[6], xaz[6], alto, azio, alto_dark, azio_dark;*/
  int32 retval, is_rise_or_set, direct;
  int32 ipl = DeterObject(ObjectName);
  helflag |= SE_HELFLAG_VISLIM_PHOTOPIC;
  /*int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);*/
  /* int32 iflag = epheflag | SEFLG_EQUATORIAL | SEFLG_TOPOCTR;*/
  if (TypeEvent == 3 || TypeEvent == 5) {
    is_rise_or_set = SE_CALC_RISE; 
    /* tret = tjdc - 3;
    if (ipl >= SE_MARS)
      tret = tjdc - 3;*/
    direct = -1;
  } else {
    is_rise_or_set = SE_CALC_SET; 
    /*tret = tjdc + 3;
    if (ipl >= SE_MARS)
      tret = tjdc + 3;*/
    direct = 1;
  }
  dtret = 999;
#if 0
  while (fabs(dtret) > 0.5) {
#else
  while (fabs(dtret) > 0.5 / 1440.0) {
#endif
    tjd += 0.7 * direct;
    if (direct < 0) tjd -= 1;
    retval = my_rise_trans(tjd, ipl, ObjectName, is_rise_or_set, helflag, dgeo, datm, &tjd, serr);
    if (retval == ERR) return ERR;
    retval = swe_vis_limit_mag(tjd, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
    if (retval == ERR) return ERR;
    while(darr[0] < darr[7]) {
      tjd += 10.0 / 1440.0 * -direct;
      retval = swe_vis_limit_mag(tjd, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
      if (retval == ERR) return ERR;
    }
    retval = time_limit_invisible(tjd, dgeo, datm, dobs, ObjectName, helflag | SE_HELFLAG_VISLIM_DARK, direct, &tret_dark, serr);
    if (retval == ERR) return ERR;
    retval = time_limit_invisible(tjd, dgeo, datm, dobs, ObjectName, helflag | SE_HELFLAG_VISLIM_NOMOON, direct, &tret, serr);
    if (retval == ERR) return ERR;
#if 0
    if (azalt_cart(tret_dark, dgeo, datm, ObjectName, helflag, darr, serr) == ERR)
      return ERR;
    if (azalt_cart(tret, dgeo, datm, ObjectName, helflag, darr+6, serr) == ERR)
      return ERR;
    dtret = acos(swi_dot_prod_unit(darr+3, darr+9)) / DEGTORAD;
#else
    dtret = fabs(tret - tret_dark);
#endif
  }
  if (azalt_cart(tret, dgeo, datm, "sun", helflag, darr, serr) == ERR)
    return ERR;
  *thel = tret;
  if (darr[1] < -12) {
    sprintf(serr, "acronychal rising/setting not available, %f", darr[1]);
    return OK;
  } else {
    sprintf(serr, "solar altitude, %f", darr[1]);
  }
  return OK;
}

static int32 get_heliacal_details(double tday, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 TypeEvent, int32 helflag, double *dret, char *serr)
{
  int32 i, retval, direct;
  AS_BOOL optimum_undefined, limit_1_undefined, limit_2_undefined;
  /* find next optimum visibility */
  optimum_undefined = FALSE;
  retval = time_optimum_visibility(tday, dgeo, datm, dobs, ObjectName, helflag, &(dret[1]), serr);
  if (retval == ERR) return ERR;
  if (retval == -2) {
    retval = OK;
    optimum_undefined = TRUE; /* change photopic <-> scotopic vision */
  }
  /* find moment of becoming visible */
  direct = 1;
  if (TypeEvent == 1 || TypeEvent == 4)
    direct = -1;
  limit_1_undefined = FALSE;
  retval = time_limit_invisible(tday, dgeo, datm, dobs, ObjectName, helflag, direct, &(dret[0]), serr);
  if (retval == ERR) return ERR;
  if (retval == -2) {
    retval = OK;
    limit_1_undefined = TRUE; /* change photopic <-> scotopic vision */
  }
  /* find moment of end of visibility */
  direct *= -1;
  limit_2_undefined = FALSE;
  retval = time_limit_invisible(dret[1], dgeo, datm, dobs, ObjectName, helflag, direct, &(dret[2]), serr);
  if (retval == ERR) return ERR;
  if (retval == -2) {
    retval = OK;
    limit_2_undefined = TRUE; /* change photopic <-> scotopic vision */
  }
  /* correct sequence of times: 
   * with event types 2 and 3 swap dret[0] and dret[2] */
  if (TypeEvent == 2 || TypeEvent == 3) {
    tday = dret[2];
    dret[2] = dret[0];
    dret[0] = tday;
    i = (int) limit_1_undefined;
    limit_1_undefined = limit_2_undefined;
    limit_2_undefined = (AS_BOOL) i;
  }
  /*if (retval == OK && dret[0] == dret[1]) */
  if (optimum_undefined || limit_1_undefined || limit_2_undefined) {
    sprintf(serr, "return values [");
    if (limit_1_undefined)
      strcat(serr, "0,");
    if (optimum_undefined)
      strcat(serr, "1,");
    if (limit_2_undefined)
      strcat(serr, "2,");
    strcat(serr, "] are uncertain due to change between photopic and scotopic vision");
  }
  return OK;
}

static int32 heliacal_ut_vis_lim(double tjd_start, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 TypeEventIn, int32 helflag, double *dret, char *serr_ret)
{
  int i;
  double d, darr[10], direct = 1, tjd, tday;
  int32 epheflag, retval = OK, helflag2;
  int32 iflag, ipl;
  int32 TypeEvent = TypeEventIn;
  char serr[AS_MAXCH];
  for (i = 0; i < 10; i++)
    dret[i] = 0;
  *dret = tjd_start;
  *serr = '\0';
  ipl = DeterObject(ObjectName);
  epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
  if (!(helflag & SE_HELFLAG_HIGH_PRECISION)) 
    iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
  if (ipl == SE_MERCURY)
    tjd = tjd_start - 30;
  else
    tjd = tjd_start - 50; /* -50 makes sure, that no event is missed, 
                         * but may return an event before start date */
  helflag2 = helflag;
  /*helflag2 &= ~SE_HELFLAG_HIGH_PRECISION;*/
  /* 
   * heliacal event
   */
  if (ipl == SE_MERCURY || ipl == SE_VENUS || TypeEvent <= 2) {
    if (ipl == -1) {
      /* find date when star rises with sun (cosmic rising) */
    retval = get_asc_obl_with_sun(tjd, ipl, ObjectName, helflag, TypeEvent, 0, dgeo, &tjd, serr);
      if (retval != OK)
	goto swe_heliacal_err; /* retval may be -2 or ERR */
    } else {
      /* find date of conjunction of object with sun */
      if ((retval = find_conjunct_sun(tjd, ipl, helflag, TypeEvent, &tjd, serr)) == ERR) {
	goto swe_heliacal_err;
      }
    }
    /* find the day and minute on which the object becomes visible */
    retval = get_heliacal_day(tjd, dgeo, datm, dobs, ObjectName, helflag2, TypeEvent, &tday, serr); 
    if (retval != OK)
      goto swe_heliacal_err;
  /* 
   * acronychal event
   */
  } else {
    if (1 || ipl == -1) {
      /*retval = get_asc_obl_acronychal(tjd, ipl, ObjectName, helflag2, TypeEvent, dgeo, &tjd, serr);*/
      retval = get_asc_obl_with_sun(tjd, ipl, ObjectName, helflag, TypeEvent, 0, dgeo, &tjd, serr);
      if (retval != OK)
	goto swe_heliacal_err;
    } else {
      /* find date of conjunction of object with sun */
      if ((retval = find_conjunct_sun(tjd, ipl, helflag, TypeEvent, &tjd, serr)) == ERR)
	goto swe_heliacal_err;
    }
    tday = tjd;
    retval = get_acronychal_day(tjd, dgeo, datm, dobs, ObjectName, helflag2, TypeEvent, &tday, serr);
    if (retval != OK)
      goto swe_heliacal_err;
  }
  dret[0] = tday;
  if (!(helflag & SE_HELFLAG_NO_DETAILS)) {
    /* more precise event times for 
     * - morning first, evening last
     * - venus and mercury's evening first and morning last
     */
    if (ipl == SE_MERCURY || ipl == SE_VENUS || TypeEvent <= 2) {
      retval = get_heliacal_details(tday, dgeo, datm, dobs, ObjectName, TypeEvent, helflag2, dret, serr);
      if (retval == ERR) goto swe_heliacal_err;
    } else if ((0)) {
      if (TypeEvent == 4 || TypeEvent == 6) direct = -1;
      for (i = 0, d = 100.0 / 86400.0; i < 3; i++, d /= 10.0) {
	while((retval = swe_vis_limit_mag(*dret + d * direct, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) == -2 || (retval >= 0 && darr[0] < darr[7])) { 
	  *dret += d * direct; 
	}
      }
      /* the last time step must be added */
      if (retval == OK)
	*dret += 1.0 / 86400.0 * direct;
    }
  } /* if (1) */
swe_heliacal_err:
  if (serr_ret != NULL && *serr != '\0')
    strcpy(serr_ret, serr);
  return retval;
}

/*###################################################################*/
static int32 moon_event_vis_lim(double tjdstart, double *dgeo, double *datm, double *dobs, int32 TypeEvent, int32 helflag, double *dret, char *serr_ret)
{
  double tjd, trise;
  char serr[AS_MAXCH];
  char ObjectName[30];
  int32 iflag, ipl, retval, helflag2, direct;
  int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
  dret[0] = tjdstart; /* will be returned in error case */
  if (TypeEvent == 1 || TypeEvent == 2) {
    if (serr_ret != NULL)
      strcpy(serr_ret, "error: the moon has no morning first or evening last");
    return ERR;
  }
  strcpy(ObjectName, "moon");
  ipl = SE_MOON;
  iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
  if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
    iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
  helflag2 = helflag;
  helflag2 &= ~SE_HELFLAG_HIGH_PRECISION;
  /* check Synodic/phase Period */
  tjd = tjdstart - 30; /* -50 makes sure, that no event is missed, 
                         * but may return an event before start date */
  if ((retval = find_conjunct_sun(tjd, ipl, helflag, TypeEvent, &tjd, serr)) == ERR)
    return ERR;
  /* find the day and minute on which the object becomes visible */
  retval = get_heliacal_day(tjd, dgeo, datm, dobs, ObjectName, helflag2, TypeEvent, &tjd, serr); 
  if (retval != OK)
    goto moon_event_err;
  dret[0] = tjd;
  /* find next optimum visibility */
  retval = time_optimum_visibility(tjd, dgeo, datm, dobs, ObjectName, helflag, &tjd, serr);
  if (retval == ERR) goto moon_event_err;
  dret[1] = tjd;
  /* find moment of becoming visible */
  /* Note: On the day of first light the moon may become visible 
   * already during day. It also may appear during day, disappear again
   * and then reappear after sunset */
  direct = 1;
  if (TypeEvent == 4)
    direct = -1;
  retval = time_limit_invisible(tjd, dgeo, datm, dobs, ObjectName, helflag, direct, &tjd, serr);
  if (retval == ERR) goto moon_event_err;
  dret[2] = tjd;
  /* find moment of end of visibility */
  direct *= -1;
  retval = time_limit_invisible(dret[1], dgeo, datm, dobs, ObjectName, helflag, direct, &tjd, serr);
  dret[0] = tjd;
  if (retval == ERR) goto moon_event_err;
#if 1
  /* if the moon is visible before sunset, we return sunset as start time */
  if (TypeEvent == 3) {
    if ((retval = my_rise_trans(tjd, SE_SUN, "", SE_CALC_SET, helflag, dgeo, datm, &trise, serr)) == ERR)
      return ERR;
    if (trise < dret[1]) {
      dret[0] = trise;
      /* do not warn, it happens too often */
      /*strcpy(serr, "start time given is sunset, but moon is observable before that");*/
    }
  /* if the moon is visible after sunrise, we return sunrise as end time */
  } else {
    if ((retval = my_rise_trans(dret[1], SE_SUN, "", SE_CALC_RISE, helflag, dgeo, datm, &trise, serr)) == ERR)
      return ERR;
    if (dret[0] > trise) {
      dret[0] = trise;
      /* do not warn, it happens too often */
      /*strcpy(serr, "end time given is sunrise, but moon is observable after that");*/
    }
  }
#endif
  /* correct order of the three times: */
  if (TypeEvent == 4) {
    tjd = dret[0];
    dret[0] = dret[2];
    dret[2] = tjd;
  }
moon_event_err:
  if (serr_ret != NULL && *serr != '\0')
    strcpy(serr_ret, serr);
  return retval;
}

static int32 MoonEventJDut(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, int32 TypeEvent, int32 helflag, double *dret, char *serr)
{
  int32 avkind = helflag & SE_HELFLAG_AVKIND;
  if (avkind) 
    return moon_event_arc_vis(JDNDaysUTStart, dgeo, datm, dobs, TypeEvent, helflag, dret, serr);
  else
    return moon_event_vis_lim(JDNDaysUTStart, dgeo, datm, dobs, TypeEvent, helflag, dret, serr);
}

static int32 heliacal_ut(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 TypeEventIn, int32 helflag, double *dret, char *serr_ret)
{
  int32 avkind = helflag & SE_HELFLAG_AVKIND;
  if (avkind) 
    return heliacal_ut_arc_vis(JDNDaysUTStart, dgeo, datm, dobs, ObjectName, TypeEventIn, helflag, dret, serr_ret);
  else
    return heliacal_ut_vis_lim(JDNDaysUTStart, dgeo, datm, dobs, ObjectName, TypeEventIn, helflag, dret, serr_ret);
}

/*' Magn [-]
' tjd_ut            start date (JD) for event search
' dgeo[3]           geogr. longitude, latitude, eye height (m above sea level)
' datm[4]           atm. pressure, temperature, RH, and VR
' - pressure        atmospheric pressure (mbar, =hPa) default 1013.25hPa
' - temperature      deg C, default 15 deg C (if at
'                   If both attemp and atpress are 0, a temperature and
'                   atmospheric pressure are estimated from the above-mentioned
'                   default values and the height above sea level.
' - RH              relative humidity in %
' - VR              VR>=1: the Meteorological range: default 40 km
'                   1>VR>0: the ktot (so the total atmospheric coefficient): 
'                   a good default would be 0.25
'                   VR=-1: the ktot is calculated from the other atmospheric 
'                   constants.
'
' dobs[6]           observer parameters
' - age [Year]      default 36, experienced sky observer in ancient times
'                   optimum age is 23
' - SN              Snellen factor of the visual aquity of the observer
'                   default 1
'                   see: http://www.i-see.org/eyecharts.html#make-your-own
' The following parameters of dobs[] are only relevant if the flag
' SE_HELFLAG_OPTICAL_PARAMS is set:
' - is_binocular    0 = monocular, 1 = binocular (actually a boolean)
' - OpticMagn       telescope magnification: 
'                   0 = default to naked eye (binocular), 1 = naked eye
' - OpticDia        optical aperture (telescope diameter) in mm
' - OpticTrans      optical transmission
'
' TypeEvent         1 morning first
'                   2 evening last
'                   3 evening first
'                   4 morning last
' dret		    output: time (tjd_ut) of heliacal event
'                   dret[0]: beginning of visibility (Julian day number)
'                   dret[1]: optimum visibility (Julian day number; 0 if SE_HELFLAG_AV)
'                   dret[2]: end of visibility (Julian day number; 0 if SE_HELFLAG_AV)
' see http://www.iol.ie/~geniet/eng/atmoastroextinction.htm
*/
int32 CALL_CONV swe_heliacal_ut(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, char *ObjectNameIn, int32 TypeEvent, int32 helflag, double *dret, char *serr_ret)
{
  int32 retval, Planet, itry;
  char ObjectName[AS_MAXCH], serr[AS_MAXCH], s[AS_MAXCH];
  double tjd0 = JDNDaysUTStart, tjd, dsynperiod, tjdmax, tadd;
  int32 MaxCountSynodicPeriod = MAX_COUNT_SYNPER;
  char *sevent[7] = {"", "morning first", "evening last", "evening first", "morning last", "acronychal rising", "acronychal setting"};
  if (dgeo[2] < SEI_ECL_GEOALT_MIN || dgeo[2] > SEI_ECL_GEOALT_MAX) {
    if (serr_ret != NULL)
      sprintf(serr_ret, "location for heliacal events must be between %.0f and %.0f m above sea\n", SEI_ECL_GEOALT_MIN, SEI_ECL_GEOALT_MAX);
    return ERR;
  }
  swi_set_tid_acc(JDNDaysUTStart, helflag, 0, serr);
  if (helflag & SE_HELFLAG_LONG_SEARCH)
    MaxCountSynodicPeriod = MAX_COUNT_SYNPER_MAX;
/*  if (helflag & SE_HELFLAG_SEARCH_1_PERIOD)
      MaxCountSynodicPeriod = 1; */
  if (serr_ret != NULL)
    *serr_ret = '\0';
  /* note, the fixed stars functions rewrite the star name. The input string 
     may be too short, so we have to make sure we have enough space */
  strcpy_VBsafe(ObjectName, ObjectNameIn);
  tolower_string_star(ObjectName);
  default_heliacal_parameters(datm, dgeo, dobs, helflag);
  swe_set_topo(dgeo[0], dgeo[1], dgeo[2]);
  Planet = DeterObject(ObjectName);
  if (Planet == SE_SUN) {
    if (serr_ret != NULL) {
      strcpy(serr_ret, "the sun has no heliacal rising or setting\n");
    }
    return ERR;
  }
  /* 
   * Moon events
   */
  if (Planet == SE_MOON) {
    if (TypeEvent == 1 || TypeEvent == 2) {
      if (serr_ret != NULL) {
        sprintf(serr_ret, "%s (event type %d) does not exist for the moon\n", sevent[TypeEvent], TypeEvent);
      }
      return ERR;
    }
    tjd = tjd0;
    retval = MoonEventJDut(tjd, dgeo, datm, dobs, TypeEvent, helflag, dret, serr);
    while (retval != -2 && *dret < tjd0) {
      tjd += 15;
      *serr = '\0';
      retval = MoonEventJDut(tjd, dgeo, datm, dobs, TypeEvent, helflag, dret, serr);
    }
    if (serr_ret != NULL && *serr != '\0')
      strcpy(serr_ret, serr);
    return retval;
  }
  /* 
   * planets and fixed stars 
   */
  if (!(helflag & SE_HELFLAG_AVKIND)) {
    if (Planet == -1 || Planet >= SE_MARS) {
      if (TypeEvent == 3 || TypeEvent == 4) {
	if (serr_ret != NULL) {
	  if (Planet == -1)
	    strcpy(s, ObjectName);
	  else
	    swe_get_planet_name(Planet, s);
	  sprintf(serr_ret, "%s (event type %d) does not exist for %s\n", sevent[TypeEvent], TypeEvent, s);
	}
	return ERR;
      }
    }
  }
  /* arcus visionis method: set the TypeEvent for acronychal events */
  if (helflag & SE_HELFLAG_AVKIND) {
    if (Planet == -1 || Planet >= SE_MARS) {
      if (TypeEvent == SE_ACRONYCHAL_RISING)
	TypeEvent = 3;
      if (TypeEvent == SE_ACRONYCHAL_SETTING)
	TypeEvent = 4;
    }
  /* acronychal rising and setting (cosmic setting) are ill-defined.
   * We do not calculate them with the "visibility limit method" */
  } else if (1) {
    if (TypeEvent == SE_ACRONYCHAL_RISING || TypeEvent == SE_ACRONYCHAL_SETTING) {
      if (serr_ret != NULL) {
	if (Planet == -1)
	  strcpy(s, ObjectName);
	else
	  swe_get_planet_name(Planet, s);
	sprintf(serr_ret, "%s (event type %d) is not provided for %s\n", sevent[TypeEvent], TypeEvent, s);
      }
      return ERR;
    }
  }
  dsynperiod = get_synodic_period(Planet);
  tjdmax = tjd0 + dsynperiod * MaxCountSynodicPeriod;
  tadd = dsynperiod * 0.6;
  if (Planet == SE_MERCURY)
    tadd = 30;
  /* 
   * this is the outer loop over n synodic periods 
   */
  tjd = tjd0;
  retval = -2;  /* indicates that another synodic period has to be done */
  for (itry = 0; 
       tjd < tjdmax && retval == -2; 
       itry++, tjd += tadd) {
    *serr = '\0';
    retval = heliacal_ut(tjd, dgeo, datm, dobs, ObjectName, TypeEvent, helflag, dret, serr);
    /* if resulting event date < start date for search (tjd0): retry starting
     * from half a period later. The event must be found now, unless there
     * is none, as is often the case with Mercury */
    while (retval != -2 && *dret < tjd0) {
      tjd += tadd;
      *serr = '\0';
      retval = heliacal_ut(tjd, dgeo, datm, dobs, ObjectName, TypeEvent, helflag, dret, serr);
    }
  }
  /* 
   * no event was found within MaxCountSynodicPeriod, return error
   */ 
  if ((helflag & SE_HELFLAG_SEARCH_1_PERIOD) && (retval == -2 || dret[0] > tjd0 + dsynperiod * 1.5)) {
    strcpy(serr, "no heliacal date found within this synodic period");
    retval = -2;
  } else if (retval == -2) {
    sprintf(serr, "no heliacal date found within %d synodic periods", MaxCountSynodicPeriod);
    retval = ERR;
  }
  if (serr_ret != NULL && *serr != '\0')
    strcpy(serr_ret, serr);
  return retval;
}