/*
 * coil.c - Create a coil of n turns between two points in space
 *
 * Alexander Enzmann
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifdef MAC
#include <console.h>
#endif
#include "def.h"
#include "lib.h"

static void
determine_coil_point(COORD4 *pos, COORD4 *norm,
		     double theta, double phi, double r0, double r1, double d0)
{
   COORD4 v0, v1, vd;
   double len;

   /* Center of the coil */
   SET_COORD4(v0, r0 * cos(theta), r0 * sin(theta), d0, 1.0);

   /* Point on the coil */
   SET_COORD4(*pos, (r0 + r1 * sin(phi)) * cos(theta),
		    (r0 + r1 * sin(phi)) * sin(theta),
		    (r1 * cos(phi)) + d0,
		    1.0);

   /* Direction from center to point */
   SUB3_COORD(vd, *pos, v0);
   len = lib_normalize_coord3(&vd);

   v0.x = r1*cos(phi)*cos(theta);
   v1.x = -(r0 + r1*sin(phi))*sin(theta);

   v0.y = r1*cos(phi)*sin(theta);
   v1.y = (r0 + r1*sin(phi))*cos(theta);

   v0.z = -r1*sin(phi);
   v1.z = 0.0;

   CROSS(*norm, v0, v1);
   len = lib_normalize_coord3(norm);

   len = DOT_PRODUCT(*norm, vd);
   if (len < 0.0) {
      norm->x *= -1.0;
      norm->y *= -1.0;
      norm->z *= -1.0;
      }
   norm->w = 0.0;
}

/* Wrap a coil around a sphere of radius sr, centered at d=t */
static void
generate_cyl_coil(COORD4 *start, COORD4 *end,
	     int turns, int step_per_turn,
	     double r0, double r1,
	     double t, double sr)
{
   int i, j, k;
   COORD4 p, p0, p1, n0;
   COORD4 dir;
   MATRIX trans, itrans;
   double dist, u0, u1, d0, d1, len;
   double deltad, deltau;
   double t0, t1, tm;
   double rt, q, rdif = r0 - r1, sr2 = sr * sr;

   SUB3_COORD(dir, *end, *start);
   dist = lib_normalize_coord3(&dir);

   /* Figure out the transform to get from the parametrically defined
      coil to the start-end line */
   lib_create_canonical_matrix(itrans, trans, start, &dir);

   /* Figure out the lower and upper limits that the sphere will
      affect the coil */
   t0 = dist * t - sr;
   t1 = dist * t + sr;
   tm = (t1 + t0) / 2.0;

   /* Calculate the step sizes for making the coil */
   deltad = dist / (double)(turns * step_per_turn);
   deltau = 2.0 * PI / (double)step_per_turn;

   /* Start generating the cylinders making up the coil */
   for (i=0,d0=0.0;i<turns;i++) {
      for (j=0;j<step_per_turn;j++,d0+=deltad) {
         d1 = d0 + deltad;
	 u0 = (double)j * deltau;
	 u1 = u0 + deltau;
	 if (d0 > t0 && d0 < t1) {
	    q = fabs(d0 - tm);
	    q = sqrt(sr2 - q * q);
	    if (q > rdif)
	       rt = r0 + (q - rdif);
	    else
	       rt = r0;
	    }
	 else
	    rt = r0;
	 determine_coil_point(&p, &n0, u0, 0.0, rt, 0.0, d0);
	 lib_transform_coord(&p0, &p, trans);
	 p0.w = r1;
	 if (d1 > t0 && d1 < t1) {
	    q = fabs(d1 - tm);
	    q = sqrt(sr2 - q * q);
	    if (q > rdif)
	       rt = r0 + (q - rdif);
	    else
	       rt = r0;
	    }
	 else
	    rt = r0;
	 determine_coil_point(&p, &n0, u1, 0.0, rt, 0.0, d1);
	 lib_transform_coord(&p1, &p, trans);
	 p1.w = r1;
	 lib_output_cylcone(&p0, &p1, OUTPUT_CURVES);
	 lib_output_sphere(&p0, OUTPUT_CURVES);
	 }
      }
   /* Cap the end */
   lib_output_sphere(&p1, OUTPUT_CURVES);
}

void
main(int argc, char *argv[])
{
    COORD4  back_color, coil_color, dir;
    COORD4  center_pt, end_pt, light;
    COORD4  from, at, up;

    MATRIX m1, m2;

   /* We are using Polyray */
   lib_set_raytracer(OUTPUT_POLYRAY);

    /* output viewpoint */
    SET_COORD(from, 0.0, 7.0,-7.0);
    SET_COORD(at,   0.0, 1.0, 2.0);
    SET_COORD(up,   0.0, 1.0, 0.0);
    lib_output_viewpoint(&from, &at, &up, 30.0, 1.0, 1.0, 512, 512);

    /* output background color - dark blue */
    SET_COORD(back_color, 0.039, 0.18, 0.376);
    lib_output_background_color(&back_color);

    /* output light source */
    SET_COORD4(light,-10.0, 10.0,-20.0, 0.7);
    lib_output_light(&light);
    SET_COORD4(light, 10.0, 10.0,-20.0, 0.7);
    lib_output_light(&light);

    /* output sphere color - blue */
    SET_COORD4(center_pt, 0.0, 1.0, 2.0, 1.5);
    SET_COORD(coil_color, 0.2, 0.2, 1.0);
    lib_output_color(NULL, &coil_color, 0.1, 0.1, 0.0, 0.4, 3.0, 0.8, 1.3);
    lib_output_sphere(&center_pt, OUTPUT_CURVES);

    /* output coil color - red */
    SET_COORD(coil_color, 1.0, 0.2, 0.2);
    lib_output_color(NULL, &coil_color, 0.2, 0.8, 0.0, 0.4, 5.0, 0.0, 0.0);

    /* compute and output coil */
    SET_COORD4(center_pt,-1.0,-1.0, 0.0, 1.0);
    SET_COORD4(end_pt, 1.0, 3.0, 4.0, 1.0);
    generate_cyl_coil(&center_pt, &end_pt, 24, 72, 1.0, 0.05, 0.5, 1.5);
}