/*
##############################################################
# COORDINATE CONVERSION ROUTINES FOR HELIOCENTRIC COORDINATES
#    part of the PDP replacement project at BBSO 1997
#
#     by: Anders Johannesson 1997-08-06
#
#    Modified:     AJO  06-Nov-97  Error Codes 
#
##############################################################
/*

/*
###########################################################
# CONVERT HELIOGRAPHIC COORDINATES INTO PROJECTED CARTESIAN
# S     - semi diameter of sun in arcseconds
# B0    - B0 in degrees
# P     - P-angle in degrees
# L0    - L0 in degrees
# B,Lin - heliographic coordinates (degrees)
# x,y   - output x,y coordinates
###################################### AJO 1997-08-06 #####
*/

#include "soleph.h"

short helio_to_cart(double S,
		    double B0,
		    double P,
		    double L0,
		    double B,
		    double Lin,
		    double *x,
		    double *y)
{

double L, Sr, B0r, Pr, L0r, Br, Lr, sina, cosrho, rho;
double sinptheta, cosptheta, ptheta, r2, r, theta;

/* CONVERT INPUT */
if (Lin < 0)
{ 
    L=360.0+Lin;
} 
else
{
    L=Lin;
}
/* L=fmod(L,360.0); */

/* CHECK INPUT */
if ( ((B >= -90) && (B <= 90))  &&  ((L >= 0) && (L <=360 )))
{

   /* CONVERSION OF PARAMETERS TO RADIANS */
   Sr=S/(3600.0*180.0/DPI);
   B0r=B0/180.0*DPI;
   Pr=P/180.0*DPI;
   L0r=L0/180.0*DPI;
   Br=B/180.0*DPI;
   Lr=L/180.0*DPI;

   /* IS THIS COORDINATE VISIBLE ON THE DISK? */
   sina=sin(Br)*sin(B0r)+cos(Br)*cos(B0r)*cos((Lin-L0)/180.0*DPI);
   if (sina <= 0)
   {
      return(FARSUN);
   }
   else
   {

      /* CALCULATIONS */
      if (fabs(B0r) > 0.0)
      {
         cosrho=(sin(Br)/cos(B0r)+cos(Br)*cos(Lr-L0r)/sin(B0r))/
    		     (1.0/tan(B0r)+tan(B0r));
      }
      else
      {
          cosrho=cos(Br)*cos(Lr-L0r);
      }
      if (cosrho >= 1.0) cosrho=1.0;
      if (cosrho <= -1.0) cosrho=-1.0;
      rho=acos(cosrho);
      if (fabs(rho) > 0.0)
      {
         sinptheta=cos(Br)*sin(Lr-L0r)/sin(rho);
         if (fabs(B0r) > 0.0)
         {
            cosptheta=(cos(rho)*cos(B0r)-cos(Br)*cos(Lr-L0r))/(sin(B0r)*sin(rho));
         }
         else
         {
            cosptheta=sin(Br)/(cos(B0r)*sin(rho));
         }
         ptheta=atan2(sinptheta,cosptheta);
      }
      else
      {
         ptheta=0.0; /* angle without meaning since r=0 */
      }

      /* RECURSIVE SOLUTION FOR r (2 ITERATIONS) */
      r2=0.0;
      r2=Sr*sin(rho+r2);
      r2=Sr*sin(rho+r2);
      r=S*sin(rho+r2);

      theta=Pr-ptheta;
      *x=-r*sin(theta);
      *y=r*cos(theta);

      return(OK);
   }
}
else
{
   if ((B < -90) || (B > 90))
      return(BADLAT);
   else
      return(BADLONG);
}

}

/*
###########################################################
# CONVERT PROJECTED CARTESIAN COORDINATES INTO HELIOGRAPHIC
# S     - semi diameter of sun in arcseconds
# B0    - B0 in degrees
# P     - P-angle in degrees
# L0    - L0 in degrees
# x,y   - location on projected disk (arcseconds)
# B,L   - output heliographic coordinates
###################################### AJO 1997-07-10 #####
*/

short cart_to_helio(double S,
		    double B0,
		    double P,
		    double L0,
		    double x,
		    double y,
		    double *B,
		    double *L)

{

double r, theta, r0, rho, sinB, sindL, cosdL, dL;
double Sr, B0r, Pr, L0r;

/* CONVERSION FROM CARTESIAN TO POLAR */
r=sqrt(x*x+y*y);
theta=atan2(-x,y);
r0=S;

/* CONVERSION OF PARAMETERS TO RADIANS */
Sr=S/(3600.0*180.0/DPI);
B0r=B0/180.0*DPI;
Pr=P/180.0*DPI;
L0r=L0/180.0*DPI;

/* CALCULATIONS */
if (r < r0)
{
   rho=asin(r/r0)-Sr*r/r0;
   sinB=cos(rho)*sin(B0r)+sin(rho)*cos(B0r)*cos(Pr-theta);
   *B=asin(sinB);
   if (fabs(fabs(*B)-DPI/2.0) > 0.0)
   {
      sindL=sin(Pr-theta)*sin(rho)/cos(*B);
      cosdL=(cos(rho)*cos(B0r)-sin(B0r)*sin(rho)*cos(Pr-theta))/cos(*B);
      dL=atan2(sindL,cosdL)+2.0*DPI;
   }
   else
   {
      dL=0.0;
   }
   *L=L0r+dL;
   *L=fmod(*L,(2.0*DPI));
   if (*L > DPI)
      *L=-2.0*DPI+*L;
   *B=*B*180.0/DPI;
   *L=*L*180.0/DPI;

   return(OK);
}
else
{
   return(BADXY);
}

}