#include "us2spherical.h"


XYZ us2rade2xyz(float ra, float dec, float radius)
/* Converts a RA (0..24) and Dec (-90..+90) to an XYZ point */
{
   XYZ xyz;
   double p, t, rsinp;

   p = (DEMAX - dec)*D2R;             	/* Convert to radians */
   t = ra*RA2R;                   	/* Convert to radians */

   rsinp =  radius * sin(p);           	/* Save one calculation */
   xyz.x =  rsinp  * cos(t);
   xyz.z = -rsinp  * sin(t);
   xyz.y =  radius * cos(p);

   return xyz;
}

void us2radept1(float ra, float dec, float radius, int sym)
{
   XYZ xyz;
   xyz = us2rade2xyz(ra, dec, radius);
   s2pt1(xyz.x, xyz.y, xyz.z, sym);
}

void us2radept(int N, float *ra, float *dec, float *radius, int *sym)
{
   int i;
   for (i=0;i<N;i++) {
      us2radept1(ra[i], dec[i], radius[i], sym[i]);
   }
}

void us2deline(float de, float ramin, float ramax, float radius)
{
   float rr = ramin*RA2R;
   float p1 = radius * sin((DEMAX-de)*D2R);
   float p2 = radius * cos((DEMAX-de)*D2R);

   float d1, d2; 
   float deg = (ramax - ramin)*15.0;
   int nseg = 100;

   XYZ centre = { 0.0, 0.0, 0.0 };
   XYZ normal = { 0.0, radius, 0.0 };
   XYZ start  = { p1, p2, 0.0 };

   /* Rotate the starting vector to ramin */
   d1 =  cos(rr) * start.x; 
   d2 = -sin(rr) * start.x; 
   start.x = d1;
   start.z = d2;

   if ((radius < 1.0) && (ramin <= 0)  && (ramax >= 24.0)) {
      s2circxz(centre.x,start.y,0, start.x, nseg, 1.0);
   } else {
      ns2varc(centre, normal, start, deg, nseg);
   }
}

void us2raline(float ra, float demin, float demax, float radius)
{
   us2radeline(ra, demin, ra, demax, radius);
}

void us2radegrid(int rastep, int destep, float radius)
{
   int i;
   for (i=RAMIN;i<=RAMAX;i+=rastep)
      us2raline((float)i, DEMIN, DEMAX, radius);

   for (i=DEMIN;i<=DEMAX;i+=destep) 
      us2deline((float)i, RAMIN, RAMAX, radius);
   
}

void us2radeline(float ra1, float de1, float ra2, float de2, float radius)
{
   XYZ xyz1, xyz2;
   XYZ centre = { 0.0, 0.0, 0.0 };
   XYZ normal;
   float angle;
   int nseg = 100;

   xyz1 = us2rade2xyz(ra1, de1, radius);
   xyz2 = us2rade2xyz(ra2, de2, radius);

   normal = cross(xyz1, xyz2);
   angle  = dota(xyz1, xyz2) * R2D;

   ns2varc(centre, normal, xyz1, angle, nseg);
} 


void us2radearrow(double ra, double de, double r1, double r2, int ci, int sz)
{  
   double p, t, rsinp; 
   float x[2], y[2], z[2];
   
   p = (90.0 - de)*D2R; 
   t = ra*15.0*D2R;
   
   rsinp = r1 * sin(p);
   x[0] = rsinp * cos(t);
   y[0] = rsinp * sin(t);
   z[0] = r1 * cos(p); 

   rsinp = r2 * sin(p);
   x[1] = rsinp * cos(t);
   y[1] = rsinp * sin(t); 
   z[1] = r2 * cos(p);
   
   s2sci(ci);
   s2sch(sz);
   s2arro(x[0],y[0],z[0],x[1],y[1],z[1]);
   s2sch(1);
}

XYZ cross(XYZ a, XYZ b)
{
   XYZ c;
   c.x = a.y*b.z - a.z*b.y;
   c.y = a.z*b.x - a.x*b.z;
   c.z = a.x*b.y - a.y*b.x;

   return c;
}

XYZ unit(XYZ a)
{
   XYZ b;
   float r;
   r = length(a);
   if (fabs(r) < SMALLEPS) {
      b.x = b.y = b.z = 0;
   } else {
      b.x = a.x/r;
      b.y = a.y/r;
      b.z = a.z/r;
   }
   return b;
}

float dotp(XYZ a, XYZ b)
{
   return (a.x*b.x + a.y*b.y + a.z*b.z);
}

float length(XYZ a)
{
   return sqrt(a.x*a.x + a.y*a.y + a.z*a.z);
}

float dota(XYZ a, XYZ b)
{
   XYZ a1 = unit(a);
   XYZ b1 = unit(b);
   float d = dotp(a1, b1);
   return acos(d);
}