/****************************************************************************** 
   This is an evolution from the first working version

   compilation:  cc pangle.c -o pangle -lm







******************************************************************************/

#define EPS 1e-3  /*  a small number  */
#define SMAX 256                                   /* max length of a string */
#define OUT stdout
#define  IN stdin
#define LOG log /* stdout */

/*  La Silla position: long =   04h 42m 55s     */
/*                      lat =  -29d 15m 24s     */
/*                     zone =  -04h         UT  */

#define PHI  -29.256667      /* Observatory latitude */


/*****************************************************************************/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

/*****************************************************************************/

/* global functions */

int Askfor(char *s, double *x, double xmin, double xmax);
void Chgest(char *ofile, char *est, char *nfile);

/*****************************************************************************/

/************************** INIZIO MAIN **************************************/

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


  FILE *log;
  int    southern =   0;              /* flag to be used later */
  char    logfile[SMAX];              /* name of log file */
  double       Pi =   3.141592653;    /* obvious */
  double     _phi = PHI;              /* Observatory latitude, temp */
  double      phi =   0;              /* Obs latitude, real */
  double        h = 100;              /* HA of the object. Set here an impossible value! */
  double      h_n = h;                /* tmp var */
  double    _delta= 100, delta=100;   /* DEC of the object. Set here an impossible value! */
  double eta, sineta, eta_, sineta_n; /* p.angle and sin of parallactic angle */
  double diffsin, eta_n, deta;        /* variables to compute    */
  double dh, deta_dh, eta_n_;         /* parallactic angle speed */

/* Create name of logfile */

  Chgest(argv[0], "log", logfile);

/*   fprintf(OUT, "\n\n Logging to: %s \n\n", logfile); */

  log = fopen(logfile, "w");

/*   If arguments given, then use them as HA and DEC */

  fprintf(LOG, "\n I have %d arguments\n", argc);

  if(argc == 3) {

    fprintf(LOG, "\n Arg 1 = %s :: Arg 2 = %s", argv[1], argv[2]);
         h = atof(argv[1]);
    _delta = atof(argv[2]);
    fprintf(LOG, "\n\n from args: Working with HA = %g and DEC = %g", h, _delta);
    fprintf(LOG,   "\n from args: And for an observatory at latitude = %g\n\n", _phi);
  }

/* Input hour angle and declination */

  Askfor("  HA ",      &h, -12.0, 12.0);
  Askfor(" DEC ", &_delta, -90.0, 90.0); 

/* Show results of read */

  fprintf(LOG, "\n\n Working with HA = %g and DEC = %g", h, _delta);
  fprintf(LOG,   "\n And for an observatory at latitude = %g\n\n", _phi);

/*   If we are in the Southern hemisphere,  */
/* work like in the Northern and change results later */

  if ( _phi < 0 ) {

    southern = 1;        _phi = -_phi;
  }

/* Convert latitude to radians		 */

  phi = _phi * Pi/180.0;

/*  --- conversion to radians */
/*  */

      h = 15 *       h * Pi / 180;
  delta =       _delta * Pi / 180;

/* --- compute sin of parallactic angle */
/* */

    sineta =  sin(h  ) *cos(phi) / sqrt(1-pow( sin(phi)*sin(delta) + cos(phi)*cos(delta)*cos(h  ),2));

/*     compute gradient in sin of parallactic angle */

  h_n = h + EPS;
  sineta_n =  sin(h_n) *cos(phi) / sqrt(1-pow( sin(phi)*sin(delta) + cos(phi)*cos(delta)*cos(h_n),2));
  diffsin = sineta_n - sineta;

/* exit if cannot compute parallactic angle */

  if(sineta > 1 || sineta < -1) { 
    fprintf(LOG, "\n\n Sorry, sin(eta)=%g, so I have to exit\n\n", sineta); 
    return -1;
  }

/*   compute p. angle, in radians, and future p. angle */

  eta = asin(sineta);
  eta_n = asin(sineta_n);

/* 		# --- conversion to degrees and hours */
/* 		# */

  eta = eta * 180 /  Pi;
  h =   h * 180   /  Pi;   /* deg   */		   
  h =   h         /  15;   /* hours */

  eta_n = eta_n * 180 /  Pi;
  h_n =   h_n * 180   /  Pi;   /* deg   */		   
  h_n =   h_n         /  15;   /* hours */

  eta_ =   eta  ; /*  for later */
  eta_n_ = eta_n; /*  for later */
		
/* --- when declination is larger than latitude, the */
/*     object passes north of the zenith, so must change */
/*     p.angle to 180 - p.angle */
/* */
/* */

  if ( _delta > _phi) { 

    eta   = h>0 ? 180 - eta   : -(180+eta  );
    eta_n = h>0 ? 180 - eta_n : -(180+eta_n);
  }


/*  --- but after p.angle reaches 90, then the value */
/*      computed by the formula is correct. So set it */
/*      back to the original value, for the corresponding */
/*      hour angle */
/*  */

    eta =  diffsin < 0 ? eta_   : eta   ;
  eta_n =  diffsin < 0 ? eta_n_ : eta_n ;
 
/* h >  $hMax ? eta_ : eta */
/* 		set eta =  h < -$hMax ? eta_ : eta */
		
/*  --- reflect around 90 deg, if Southern hemisphere */
/*  */

  if ( southern == 1 ) {

    eta   = h > 0 ? 180 - eta   : -eta  -180 ;
    eta_n = h > 0 ? 180 - eta_n : -eta_n-180 ;
  }

/* compute rotation speed of parallactic angle */

  deta = eta_n - eta;
  dh = h_n - h;
  deta_dh = deta/dh;

  fprintf(LOG, "\n At HA = %g I have sin(eta) = %g => eta = %g diff = %g :: new eta and new h are %g and %g so speed is %g deg/hr\n\n", h, sineta, eta, diffsin, eta_n, h_n, deta_dh);
  fprintf(OUT, "%g %g %g %g\n", h, _delta, eta, deta_dh);

  fclose(log);

  return 1;
}

/* ========================================================================= */
/* ========================================================================= */
/* ========================================================================= */
/* ========================================================================= */

/* Function to read variables with constrains on their value */

  int Askfor(char *s, double *x, double xmin, double xmax) {

    while(*x>xmax || *x<xmin) {

      fprintf(OUT, "\n Remember you can call me as 'pangle HA DEC'\n\n");
      fprintf(OUT, " Give %s between %g and %g : ", s, xmin, xmax);
      fscanf(IN, "%lf", x);
    }
    return 0;
  }


void Chgest(char *ofile, char *est, char *nfile) {

	int i = 0, j = 0, punto = 0, c;

/* --- Anzitutto controlla se nel vecchio nome c'e' il punto
*/
	while((c = ofile[i]) != '\0') {

		if(c == '.') {

			punto++;
		}
		i++;
	}
	i = 0;

/* --- Se c'e' il punto assegna la prima parte del nome del file
*/
	if(punto) {

		while(ofile[i] != '.') {

			nfile[i] = ofile[i];
			i++;
		} 
	} else {

		while(ofile[i] != '\0') {

			nfile[i] = ofile[i];
			i++;
		} 
	}
	nfile[i] = '.';
	i++;

	while(est[j] != '\0') {

		nfile[i] = est[j];
		i++;
		j++;
	}
	nfile[i] = '\0';
}

/**************************************************************************
	Dato in input un nome di file da' in output una
        stringa contenente il nome del file senza estensione     */

void ExtEst(char *ofile, char *nfile) {

	int i = 0, j = 0, punto = 0, c;

/* --- Anzitutto controlla se nel vecchio nome c'e' il punto
*/
	while((c = ofile[i]) != '\0') {

		if(c == '.') {

			punto++;
		}
		i++;
	}
	i = 0;

/* --- Se c'e' il punto assegna la prima parte del nome del file
*/
	if(punto) {

		while(ofile[i] != '.') {

			nfile[i] = ofile[i];
			i++;
		} 
	} else {

		while(ofile[i] != '\0') {

			nfile[i] = ofile[i];
			i++;
		} 
	}
	nfile[i] = '\0';
}

/* ========================================================================= */
/* ========================================================================= */
/* ========================================================================= */
/* ========================================================================= */