//  three.c provides contour plots in projections, c-splines, and some miscellaneous stuff.

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


//  Contour plots

  static void contourline(contourtask *CT, const vector *V) {
if (CT->pen) vectlineto(CT->wind,CT->proj,V); else vectmoveto(CT->wind,CT->proj,V); CT->pen=1; }

  static bool advancecontour(contourtask *CT, void **L, void **R) {
void *N; vector V; if (!L || !*L || !R || !*R) return 0;
N=CT->neigh(*L,*R); if (!N) return 0;
CT->where(N,&V); if (V.z<CT->level) *R=N; else *L=N; return 1; }

  static void segmoveline(contourtask *CT, void *L, void *R) {
vector VR,VL,V; real Frac; if (!L || !R) return;
CT->where(R,&VR); CT->where(L,&VL);
Frac= (VL.z - CT->level) / (VL.z - VR.z);
V.x= (Frac * VR.x) + ((1.0 - Frac) * VL.x);
V.y= (Frac * VR.y) + ((1.0 - Frac) * VL.y);
V.z=CT->level; contourline(CT,&V); }

  void drawcontour(contourtask *CT, real Level) {
int count; void *L,*R; char *S; vector V; image D; window *W;
if (!CT) return; W=CT->wind; if (!W) return;
CT->level=Level; L=CT->start(CT->level); if (!L) return; R=CT->loop(L);
CT->pen=0; pencolour(W,&(CT->rgb)); beginpath(W);
for (count=0;count<CT->max;count++) {
  segmoveline(CT,L,R);
  if (CT->done(L,R) || !advancecontour(CT,&L,&R)) break; }
stroke(W); penindex(W,black); if (!CT->ticks) return;
count*=CT->tickfrac;
L=CT->start(CT->level); if (!L) return; R=CT->loop(L);
fillcolour(W,&(CT->rgb));
for (;count>=0;count--) { if (CT->done(L,R) || !advancecontour(CT,&L,&R)) break; }
//CT->pen=0; segmoveline(CT,L,R);
CT->where(R,&V); projectpoint(CT->proj,&V,&D);
S=fixedstring(Level,3,2); textcentre(W,D.pt.x,D.pt.y - 5,S); old(S); fillindex(W,black); }

  void framecontourplot(contourtask *CT) {
int count; void *P; vector V; window *W;
if (!CT) return; W=CT->wind; if (!W) return;
P=CT->outline; CT->pen=0; pencolour(W,&(CT->rgb)); beginpath(W);
for (count=0;count<CT->max;count++) { CT->where(P,&V);
  contourline(CT,&V); P=CT->loop(P); if (P==CT->outline) break; }
CT->where(P,&V); contourline(CT,&V); stroke(W); penindex(W,black); }


//  N-by-N matrices

//  This procedure solves an N by N tridiagonal matrix.  A[] are the elements to the
//  left of the diagonal, B[] are the diagonal elements, and C[] are to the right of
//  it.  A[0] and C[N - 1] are ignored.  A[], B[] and C[] contain garbage on exit.
//  D[] are the driving terms on entry, and contain the solution on exit.

  void solvetridiag(int N, real *A, real *B, real *C, real *D) {
int i; real mult; for (i=1;i<N;i++) { mult= A[i] / B[i - 1];
B[i]-= C[i - 1] * mult; D[i]-= D[i - 1] * mult; }
for (i= N - 2;i>=0;i--) D[i]-= D[i + 1] * C[i] / B[i + 1];
for (i=0;i<N;i++) D[i]/=B[i]; }


//  C-splines

/*  void drawcspline(int OFFSET, real *H, real *V, real *HT, real *VT, int SEGS) {
int i; real t,t2,B2H,B3H,B2V,B3V;
H+=OFFSET; V+=OFFSET; HT+=OFFSET; VT+=OFFSET;
B2H= (3.0 * (H[1] - H[0])) - (2.0 * HT[0]) - HT[1];
B3H= (2.0 * (H[0] - H[1])) +        HT[0]  + HT[1];
B2V= (3.0 * (V[1] - V[0])) - (2.0 * VT[0]) - VT[1];
B3V= (2.0 * (V[0] - V[1])) +        VT[0]  + VT[1];
//if (OFFSET==0) { colour(pink); drawcircle(H[0],V[0],2); // LineTo(H[1],V[1]);
// colour(grey); }
MoveTo(H[0],V[0]); for (i=1;i<SEGS;i++) { int h,v; t= i / (real)SEGS; t2= t * t;
h= H[0] + (HT[0] * t) + (B2H * t2) + (B3H * t * t2);
v= V[0] + (VT[0] * t) + (B2V * t2) + (B3V * t * t2);
LineTo(h,v); } LineTo(H[1],V[1]); }*/


//  Platonic solids

const static char DodA[]="1357911233345556777899901224446668880002";
const static char DodB[]="1111112321232123212321232344344344344344";

  void framedodecahedron(window *W, vector *Offset) {
projection *Proj; real Phi,Theta; vector C,V; int I; image D; bool Pen=0;
if (!W) return; Proj=W->proj; if (!Proj) return;
if (Offset) C=*Offset; else zerovector(&C); beginpath(W);
for (I=0;I<40;I++) {
  Phi= ((real)(DodA[I] - '0')) * pi / 5.0; switch (DodB[I]) {
    case '1': Theta=0.6523581; break;
    case '2': Theta=1.3820859; break;
    case '3': Theta=1.7595068; break;
    case '4': Theta=2.4892346; break; }
  polartocartesian(1.0,Theta,Phi,&V); vectorincr(&V,&C);
  projectpoint(Proj,&V,&D); movetolineto(W,&D,&Pen); }
stroke(W); }

const static char BCCA[]="100332330012";
const static char BCCB[]="321123455433";

  void framebccwignerseitzcell(window *W, vector *Offset) {
projection *Proj; vector V; int I,J; image D; bool Pen;
if (!W) return; Proj=W->proj; if (!Proj) return; beginpath(W);
for (J=0;J<4;J++) { Pen=0; for (I=0;I<12;I++) {
  switch (BCCA[I]) {
    case '0': V.x=1.0; V.y=0.0; break;
    case '1': V.x=1.0; V.y=0.5; break;
    case '2': V.x=0.5; V.y=1.0; break;
    case '3': V.x=0.0; V.y=1.0; break; }
  switch (BCCB[I]) {
    case '1': V.z= 1.0; V.x/=2.0; V.y/=2.0; break;
    case '2': V.z= 0.5; break;
    case '3': V.z= 0.0; break;
    case '4': V.z=-0.5; break;
    case '5': V.z=-1.0; V.x/=2.0; V.y/=2.0; break; }
  switch (J) { case 1: rotz090(&V); break; case 2: rotz180(&V); break; case 3: rotz270(&V); }
  if (Offset) vectorincr(&V,Offset);
  projectpoint(Proj,&V,&D); movetolineto(W,&D,&Pen); } }
stroke(W); }

const static char FCCA[]="212123232";
const static char FCCB[]="123454321";

  void framefccwignerseitzcell(window *W, vector *Offset) {
projection *Proj; vector V; int I,J; image D; bool Pen;
if (!W) return; Proj=W->proj; if (!Proj) return; beginpath(W);
for (J=0;J<4;J++) { Pen=0; for (I=0;I<9;I++) {
  switch (FCCA[I]) {
    case '1': V.x=0.5; V.y=-0.5; break;
    case '2': V.x=1.0; V.y= 0.0; break;
    case '3': V.x=0.5; V.y= 0.5; break; }
  switch (FCCB[I]) {
    case '1': V.z= 1.0; V.x=0.0; V.y=0.0; break;
    case '2': V.z= 0.5; break;
    case '3': V.z= 0.0; break;
    case '4': V.z=-0.5; break;
    case '5': V.z=-1.0; V.x=0.0; V.y=0.0; break; }
  switch (J) { case 1: rotz090(&V); break; case 2: rotz180(&V); break; case 3: rotz270(&V); }
  if (Offset) vectorincr(&V,Offset);
  projectpoint(Proj,&V,&D); movetolineto(W,&D,&Pen); } }
stroke(W); }