//  sphcoord.c
//
//  Thuis is a generalised group of routines for drawing coordinate line on
//  a spherical surface, be they lines of latitude & longitude, or RA & Dec.
//  It can be used for perspective views, or for the surface of spherical objects.
//
//  by John D. de Boer

#include "common.h"
#include "cstring.h"
#include "cyber.h"
#include "disk.h"
#include "floating.h"
#include "frame.h"
#include "graph.h"
#include "gui.h"
#include "vector.h"
#include "view.h"
#include "sphcoord.h"


const real MinsInRad = 3437.7467708;  //  the number of minutes of arc in one radian


//  Bounds of perspective views

  static void station(const projection *Proj, int i, point *P) {
if (!Proj || !P) return;
P->x= left(Proj->port) + ((i % 3) * 0.5 * width(Proj->port));
P->y= top(Proj->port) - ((i / 3) * 0.5 * height(Proj->port)); }

  static real lsr3(real *X, int A, int B, int C) {
real R; R=X[A]; if (X[B]<R) R=X[B]; if (X[C]<R) R=X[C]; return R; }

  static real gtr3(real *X, int A, int B, int C) {
real R; R=X[A]; if (X[B]>R) R=X[B]; if (X[C]>R) R=X[C]; return R; }


//  This function determines what portion of a sphere is visible in the view
//  rectangle.  The maximum bounds are the coordinate lines which have any part
//  visible.  The minimum bounds are the coordinate lines which cross the view
//  completely, vertically or horizontally, or which don't go outside the view
//  rectangle.

  static void boundsofview(const projection *Proj, bool Globe, real Margin, spherebounds *B) {
int i; point P; vector V; image D;
real MidPhi,    //  the longitude at the centre of the view
     Theta[9],  //  the polar coordinates corresponding to 9 points in grid on the screen.
     Phi[9];    //    Point #0 is at the top left, and they go in rows, starting across.
bool LtoR,      //  whether azimuth increases going from left to right across screen
     Out,       //  whether any point falls off the sphere
     Good;      //  whether a reverse projection succeeds

zerovector(&V);
V.z= 1.0; if (Globe) projontosphere(Proj,&V,&D); else projectdirection(Proj,&V,&D);
B->north=(D.good && D.miss<=0);
V.z=-1.0; if (Globe) projontosphere(Proj,&V,&D); else projectdirection(Proj,&V,&D);
B->south=(D.good && D.miss<=0);

LtoR=Globe; Out=0; MidPhi=modtwopisym(Proj->azim + (LtoR ? pi : 0.0));

for (i=0;i<9;i++) {
  station(Proj,i,&P);
  Good= Globe ? reverseprojontosphere(Proj,P,&V) : reverseprojection(Proj,P,&V);
  if (!Good) { Theta[i]= pibytwo - Proj->elev; Phi[i]=MidPhi; Out=1; continue; }
  cartesiantopolar(&V,NULL,Theta + i,Phi + i); Phi[i]=modtwopisym(Phi[i]);
  if (i % 3 == (LtoR ? 0 : 2) && Phi[i]>MidPhi) Phi[i]-=twopi;
  if (i % 3 == (LtoR ? 2 : 0) && Phi[i]<MidPhi) Phi[i]+=twopi; }

B->max.north=lsr3(Theta,0,1,2);
B->max.south=gtr3(Theta,6,7,8);
B->max.west= LtoR ? lsr3(Phi,0,3,6) : lsr3(Phi,2,5,8);
B->max.east= LtoR ? gtr3(Phi,2,5,8) : gtr3(Phi,0,3,6);

for (i=0;i<9;i++) {
  station(Proj,i,&P);
  if (i<=2) P.y-=Margin; if (i>=6) P.y+=Margin;
  if (i % 3 == 0) P.x+=Margin; if (i % 3 == 2) P.x-=Margin;
  Good= Globe ? reverseprojontosphere(Proj,P,&V) : reverseprojection(Proj,P,&V);
  if (!Good) { Theta[i]= pibytwo - Proj->elev; Phi[i]=MidPhi; Out=1; continue; }
  cartesiantopolar(&V,NULL,Theta + i,Phi + i); Phi[i]=modtwopisym(Phi[i]);
  if (i % 3 == (LtoR ? 0 : 2) && Phi[i]>MidPhi) Phi[i]-=twopi;
  if (i % 3 == (LtoR ? 2 : 0) && Phi[i]<MidPhi) Phi[i]+=twopi; }

B->min.north=gtr3(Theta,0,1,2);
B->min.south=lsr3(Theta,6,7,8);
B->min.west= LtoR ? gtr3(Phi,0,3,6) : gtr3(Phi,2,5,8);
B->min.east= LtoR ? lsr3(Phi,2,5,8) : lsr3(Phi,0,3,6); 
//  The coord line tangent to the edge of the screen can actually occur at any point along the side;
//  3 points is a kludge.

if (B->north || Out) B->max.north=B->min.north=0.0;
if (B->south || Out) B->max.south=B->min.south=pi;
if (B->north || B->south || Out) { B->max.west=-pi; B->max.east=pi; }

 }

/*  static void testbounds(const projection *Proj, spherebounds *B) {
colour(green);
MoveTo(100,100); drawstr("N:   "); drawint(B->max.north * fiftyseven); drawstr("¡");
MoveTo(100,115); drawstr("S:   "); drawint(B->max.south * fiftyseven); drawstr("¡");
MoveTo(100,135); drawstr("W:   "); drawint(B->max.west * fiftyseven); drawstr("¡");
MoveTo(100,150); drawstr("mid: "); drawint(modtwopisym(Proj->azim + pi) * fiftyseven); drawstr("¡");
MoveTo(100,165); drawstr("E:   "); drawint(B->max.east * fiftyseven); drawstr("¡"); }*/


//  Trig look-up table

bool TrigInitted = 0;  //
real *SinLkUp = NULL,    //  sine and cosine look-up table; initialised by initspheretask();
     *CosLkUp = NULL;

  static void inittriglookup(void) {
int i; SinLkUp=(real *)pointer(360 * sizeof(real));
for (i=0;i<360;i++) SinLkUp[i]=sin(i / fiftyseven);
CosLkUp=pointer(360 * sizeof(real));
for (i=0;i<360;i++) CosLkUp[i]=cos(i / fiftyseven);
TrigInitted=1; }

static real sinlu(long DEG) { if (DEG<0) DEG+=360; return SinLkUp[DEG]; }

static real coslu(long DEG) { if (DEG<0) DEG+=360; return CosLkUp[DEG]; }


//  This routine sets default values for all of the "spheretask" parameters.

  void setspheretaskradius(spheretask *ST, real F, bool Hour) {  //  F is the radius in pixels
/**/ if (F>5000.0) { ST->dlat=60; ST->dlong= Hour ? 75 : 60; ST->mlat= ST->mlong= 300; }
else if (F>2500.0) { ST->dlat=120; ST->dlong= Hour ? 150 : 120; ST->mlat= ST->mlong=600; }
else if (F>1000.0) { ST->dlat= ST->dlong= 300; ST->mlat=600; ST->mlong=900; }
else if (F> 200.0) { ST->dlat=600; ST->dlong=900; ST->mlat= ST->mlong= 1800; }
else               { ST->dlat= ST->mlat= ST->dlong= ST->mlong= 1800; }
ST->latintv=ST->dlat; ST->longintv=ST->dlong; }

  void initspheretask(spheretask *ST, window *W, projection *P, bool Globe, real Margin, sphmapptr Map, bool Hour) {
real F,TickFactor; if (!ST || !W || !P) return;
if (!TrigInitted) inittriglookup();
ST->wind=W; ST->proj=P; boundsofview(P,Globe,Margin,&(ST->bound));  //  args should be 1,5,...,0 for Twistory
ST->map=Map; ST->hour=Hour;
setcolour(&(ST->c1st),white); setcolour(&(ST->c2nd),white);
ST->lat1= 5400 - (ST->bound.max.south * MinsInRad);
ST->lat2= 5400 - (ST->bound.max.north * MinsInRad);
ST->long1= ST->bound.max.west * MinsInRad;
ST->long2= ST->bound.max.east * MinsInRad;

F= P->persp ? P->pix : P->mag;

setspheretaskradius(ST,F,Hour);

/*     if (F>5000.0) { ST->dlat=60; ST->dlong= Hour ? 75 : 60; ST->mlat= ST->mlong= 300; }
else if (F>2500.0) { ST->dlat=120; ST->dlong= Hour ? 150 : 120; ST->mlat= ST->mlong=600; }
else if (F>1000.0) { ST->dlat= ST->dlong= 300; ST->mlat=600; ST->mlong=900; }
else if (F> 200.0) { ST->dlat=600; ST->dlong=900; ST->mlat= ST->mlong= 1800; }
else               { ST->dlat= ST->mlat= ST->dlong= ST->mlong= 1800; }

ST->latintv=ST->dlat; ST->longintv=ST->dlong;*/

TickFactor= 1.7 * cos(P->elev); if (TickFactor<1.0) {
       if (TickFactor * F > 5000.0) ST->longintv= Hour ?  75 :  60;
  else if (TickFactor * F > 2500.0) ST->longintv= Hour ? 150 : 120;
  else if (TickFactor * F > 1000.0) ST->longintv= 300;
  else if (TickFactor * F >  500.0) ST->longintv= Hour ? 450 : 600;
  else if (TickFactor * F >  200.0) ST->longintv= 900;
  else                              ST->longintv=1800; }

if (ST->bound.min.south<pibytwo) ST->ticklat= 5400.0 - (MinsInRad * ST->bound.min.south);
else if (ST->bound.min.north>pibytwo) ST->ticklat= 5400.0 - (MinsInRad * ST->bound.min.north);
else ST->ticklat=0.0;
ST->ticklong= MinsInRad * (P->persp ? ST->bound.min.east : ST->bound.min.west); }


//  Utility routines

  static long longitudeshortening(spheretask *ST, long LONG) {
if (LONG %    5) return 300;  //  if not an even  5' or 30 sec
if (LONG %   15) return 150;  //  if not an even 15' or a whole min
if (ST->hour) {
  if (LONG % 75) return  30;  //  if not an even 5 min
} else {
  if (LONG % 30) return  60;  //  if not an even 30'
  if (LONG % 60) return  30;  //  if not a whole degree
}
if (LONG %  300) return  10;  //  if not an even  5¡ or 20 min
if (LONG %  900) return   4;  //  if not an even 15¡ or a whole hour
if (LONG % 1800) return   2;  //  if not an even 30¡ or 2 hours
if (LONG % 5400) return   1;  //  if not an even 90¡ or 6 hours
return 0; }  //  every 90¡, the meridian goes right to the pole

  static long roundstartbound(long START, long BOUND, long FACTOR) {
long i; i=gtr(START,BOUND); return i - mod(i,FACTOR); }

  static long roundendbound(long END, long BOUND, long FACTOR) {
long i,j; i=lsr(END,BOUND); j=mod(i,FACTOR); if (j) i+= FACTOR - j; return i; }


//

  static int circleresolution(int SPACING) {  //  argument in minutes; return in degrees
if (SPACING > 5 * 60) return 5;
if (SPACING > 2 * 60) return 3;
return SPACING / 60; }

  static void drawlat(spheretask *ST, long LAT) {  //  can be sped up.  draws entire parallels
long LONG,STEP; real Theta,SinT,CosT; vector U; image D; bool PenDown; window *W;
W=ST->wind; if (!W) return;
PenDown=0; Theta= (5400 - LAT) / MinsInRad; SinT=sin(Theta); CosT=cos(Theta);
STEP=circleresolution(ST->dlong);
beginpath(W);
for (LONG=-180;LONG<=180;LONG+=STEP) {
  U.x= SinT * coslu(LONG); U.y= SinT * sinlu(LONG); U.z=CosT;
  ST->map(ST->proj,&U,&D); movetolineto(W,&D,&PenDown); }
stroke(W); }

  static void drawlong(spheretask *ST, long LONG, long LAT1, long LAT2) {
long THETA,SKIP,STEP,T1,T2; real Phi,SinP,CosP,SinT; vector U; image D; bool PenDown,Last; window *W;
W=ST->wind; if (!W) return;
SKIP= ST->dlat * longitudeshortening(ST,LONG); SKIP/=60;  //  might not work if dlat<1¡
PenDown=0; Phi= LONG / MinsInRad; SinP=sin(Phi); CosP=cos(Phi);
STEP=circleresolution(ST->dlat);
T1= (5400 - LAT2) / 60; T2= (5400 - LAT1) / 60;
beginpath(W);
for (THETA=SKIP;THETA <= 180 + STEP;THETA+=STEP) {
  if (THETA + (2 * STEP) < T1) continue;
  if (THETA < 180 - SKIP) Last=0; else { Last=1; THETA= 180 - SKIP; }
  SinT=sinlu(THETA); U.x= SinT * CosP; U.y= SinT * SinP; U.z=coslu(THETA);
  ST->map(ST->proj,&U,&D); movetolineto(ST->wind,&D,&PenDown);
  if (Last || THETA>T2) { stroke(W); return; } }
stroke(W); }

  void drawsphericalcoords(spheretask *ST) {
long LAT1,LAT2,    //  rounded latitude range
     LONG1,LONG2;  //  rounded longitude range
long i; window *W; if (!ST || !ST->map) return;
W=ST->wind; if (!W) return;

if (ST->dlat<ST->mlat) { pencolour(W,&(ST->c2nd));
LAT1=roundstartbound(ST->lat1,-5400,ST->dlat);
LAT2=roundendbound(ST->lat2,5399,ST->dlat);
LONG1=roundstartbound(ST->long1,-21600 + ST->dlong,ST->dlong);
LONG2=roundendbound(ST->long2,21600,ST->dlong);
for (i=LAT1; i<=LAT2 ;i+=ST->dlat ) if (i % ST->mlat) drawlat(ST,i);
for (i=LONG1;i<=LONG2;i+=ST->dlong) if (i % ST->mlong) drawlong(ST,i,LAT1,LAT2); }

pencolour(W,&(ST->c1st));
LAT1=roundstartbound(ST->lat1,-5400,ST->mlat);
LAT2=roundendbound(ST->lat2,5399,ST->mlat);
LONG1=roundstartbound(ST->long1,-21600 + ST->mlong,ST->mlong);
LONG2=roundendbound(ST->long2,21600,ST->mlong);
for (i=LAT1; i<=LAT2 ;i+=ST->mlat ) drawlat(ST,i);
for (i=LONG1;i<=LONG2;i+=ST->mlong) drawlong(ST,i,LAT1,LAT2);
penindex(W,black); }


//  Drawing special circles

  void drawequator(spheretask *ST) {
if (!ST || !ST->map) return; drawlat(ST,0); }

  void drawdashedequator(spheretask *ST) {
window *W; projection *Proj; float Len,Dash[2];
long CENTRE,LONG,STEP; vector U; image D; bool PenDown;
if (!ST || !ST->map) return;
W=ST->wind; if (!W) return; Proj=W->proj; if (!Proj) return;

Len=6.0; Dash[0]=Len; Dash[1]=Len;
CENTRE= fiftyseven * Proj->azim; CENTRE-= CENTRE % 15;
PenDown=0;
STEP=circleresolution(ST->dlong);
savegraphics(W); CGContextSetLineDash(W->q2d,Len / 2,Dash,2);
beginpath(W);
for (LONG=CENTRE;LONG <= CENTRE + 180;LONG+=STEP) {
  U.x=cos(LONG / fiftyseven); U.y=sin(LONG / fiftyseven); U.z=0.0;
  ST->map(Proj,&U,&D); movetolineto(W,&D,&PenDown); }
stroke(W); beginpath(W);
for (LONG=CENTRE;LONG >= CENTRE - 180;LONG-=STEP) {
  U.x=cos(LONG / fiftyseven); U.y=sin(LONG / fiftyseven); U.z=0.0;
  ST->map(Proj,&U,&D); movetolineto(W,&D,&PenDown); }
stroke(W);
restoregraphics(W); }

  void drawprimemeridian(spheretask *ST, bool Colure) {
if (!ST || !ST->map) return;
drawlong(ST,0,-5400,5400); if (Colure) drawlong(ST,10800,-5400,5400); }


//  Determining locations for tick marks

  void coordtick(spheretask *ST, real LAT, real LONG, image *D, char Dir, vector *U) {
vector V,W; image E; char C[2]; if (!ST || !ST->map || !D) return;
polartocartesian(1.0,(5400 - LAT) / MinsInRad,LONG / MinsInRad,&V); ST->map(ST->proj,&V,D);
if (!U) return; C[0]=Dir; C[1]=0; makelowercase(C);
switch (C[0]) {
  case 'n': LAT+=ST->dlat; break;
  case 'e': LONG+=ST->dlong; break;
  case 'w': LONG-=ST->dlong; break;
  case 's': default: LAT-=ST->dlat; }
polartocartesian(1.0,(5400 - LAT) / MinsInRad,LONG / MinsInRad,&W); ST->map(ST->proj,&W,&E);
if (!D->good || !E.good) { zerovector(U); return; }
setvector(U,E.pt.x - D->pt.x,E.pt.y - D->pt.y,0.0); normalise(U); }