#include "elastro.h"





/*
  -----------------------------------------------------------------------
  This sourcecode is stolen from Elwood Downey's gratefully
  program-sources XEphem Vers. 3.5.2.
  It includes the routines from XEphem refract.c

  Because i'm not the mathematical hero, i use it for my own
  program.
  I've asked him, if this is possible, and i'm very happy,
  that he says yes.
 
  -------------------------------------------------------- 
  Many thanks to Elwood, - also for his wonderfull XEphem.
  --------------------------------------------------------
  
  But, have a look, this source are under the same licence-agreements,
  as all the others.
  You can use it by your own, but you cannot distribute it, nor
  go any other way to earn money with this.
  Because i love the program XEphem, and my program will use it's
  features, i use the same licence-agreements.

  The modifications i've done, are only for code-design, because
  Elwood are one of the great-old-man programmers, which are
  100 percent Kernigan-Richie compatible.
  I hope he is not angry, about my changes of design, because the
  only reason for that is my own, subjective, aversion against
  style-guides.

  If you will see proper sources, get the sources for XEphem,
  (install, and use it :-), and have a look for the sourcefiles
  you found there.
  -------------------------------------------------------------------------
/*






/* ###############################################
   REFRACTION #################################### */





static void unrefract ( double pr, double tr, double aa, double * ta)
{

#define	LTLIM	14.5
#define	GELIM	15.5

	double aadeg = deg(aa);

	if (aadeg < LTLIM) {
	    unrefractLT15 (pr, tr, aa, ta);
	}
	else if (aadeg >= GELIM) {
	    unrefractGE15 (pr, tr, aa, ta);
	}
	else {
	    /* smooth blend -- important for inverse */
	    double taLT, taGE, p;

	    unrefractLT15 (pr, tr, aa, &taLT);
	    unrefractGE15 (pr, tr, aa, &taGE);
	    p = (aadeg - LTLIM)/(GELIM - LTLIM);
	    *ta = taLT + (taGE - taLT)*p;
	}
}







static void unrefractGE15 ( double pr, double tr, double aa, double * ta)
{
	double r;
	
	r = 7.888888e-5*pr/((273+tr)*tan(aa));
	*ta  =  aa - r;
}








static void unrefractLT15 ( double pr, double tr, double aa, double * ta)
{
	double aadeg = deg(aa);
	double r, a, b;

	a = ((2e-5*aadeg+1.96e-2)*aadeg+1.594e-1)*pr;
	b = (273+tr)*((8.45e-2*aadeg+5.05e-1)*aadeg+1);
	r = rad(a/b);

	*ta  =  (aa < 0 && r < 0) ? aa : aa - r;	/* 0 below ~5 degs */
}









/* correct the true altitude, ta, for refraction to the apparent altitude, aa,
 * each in radians, given the local atmospheric pressure, pr, in mbars, and
 * the temperature, tr, in degrees C.
 */
static void refract ( double pr, double tr, double ta, double * aa)
{
#define	MAXRERR	rad(0.1/3600.)	/* desired accuracy, rads */

	double d, t, t0, a;

	/* first guess of error is to go backwards.
	 * make use that we know delta-apparent is always < delta-true.
	 */
	unrefract (pr, tr, ta, &t);
	d = 0.8*(ta - t);
	t0 = t;
	a = ta;

	/* use secant method to discover a value that unrefracts to ta.
	 * max=7 ave=2.4 loops in hundreds of test cases.
	 */
	while (1) {
	    a += d;
	    unrefract (pr, tr, a, &t);
	    if (fabs(ta-t) <= MAXRERR)
		break;
	    d *= -(ta - t)/(t0 - t);
	    t0 = t;
	}

	*aa = a;

#undef	MAXRERR
}