///////////////////////////////////////////////////////////////
//
// K3-000d.cpp
//
//   arc patch in ph=pi/2, rho=0
//
// notation: ijk=phase label, 64 cases i=0:3
//   a,b,c,d,e,f,g,h,i,j... = skeleton surface
//    patches of 4-polytope for one phase.
//
// - Define the equation of a quartic  Calabi-Yau K3
//   A. Hanson Oct 2005   
//
// original version by Philip Fu (cwfu@acm.org)
//
// 3/9/2003 1:18PM
//
// All rights reserved
//
// NOTE:
// - To add an equation, make modifications marked by "EDIT HERE"
//
///////////////////////////////////////////////////////////////


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "ndeqdef.h"
#include "eqcommon.h"


#include "K3.h"   // Parameters XIMAX, XISTEPS, XISCALE, THSTEPS

///////////////////////////////////////////////////////////////
// (PART A) MATHEMATICAL NOTE
///////////////////////////////////////////////////////////////

// *** <0> EDIT HERE (optional)

#define K3_000_DESC "\
\n\
  <Quartic CP3 K3 - 000>\n\
\n\
  - 2 parameters: theta1 and theta2 range from 0 to 2Pi\n\
  -  1 complex cubic equation:\n\
  - It is a 2D manifold embedding in 4D space\n"

// POWERS of general equation can be substituted HERE

#define Nval  (4.0)
#define K1val  (0.0)
#define K2val  (0.0)
#define K3val  (0.0)


///////////////////////////////////////////////////////////////
// (PART B) INDIVIDUAL FUNCTIONS FOR EACH DIMENSION
///////////////////////////////////////////////////////////////



////////////////////////////////////////////////////
// Number of Parameters and Dimensions
////////////////////////////////////////////////////

// *** <1> EDIT HERE

// number of parameters
#define NUM_PARA	2

// number of dimensions (embedding dimension)
#define NUM_DIM		16



////////////////////////////////////////////////////
// You may give name to the parameters by #define
////////////////////////////////////////////////////

// *** <2> EDIT HERE (optional)

//#define th		(in_para_array[0])
//#define xi		(in_para_array[1])

 #define ph		(in_para_array[0])
 #define rh		(in_para_array[1])

#define DIM(i)		(out_value_array[(i)])

#define th  (0.0)
#define xi (0.0)


////////////////////////////////////////////////////
// Define the actual equations defined here
////////////////////////////////////////////////////

// NOte:
// - NUM_PARA - number of elements in in_para_array
//   NUM_DIM  - number of elements in out_value_array
// - Make sure the array index isn't out of range, otherwise, may crash!

// *** <3> EDIT HERE


PLUGIN_API void
computeEquation( const double in_para_array [_MAX_PARAM] , double out_value_array [_MAX_DIM] )
{

  double xx,yy,zz,ww,rr,rt2,X0,Y0,X1,Y1,X2,Y2,X3,Y3,denom,expxi,exprh,angle;
  rt2 = sqrt(2.0);
  expxi = XISCALE*sinh(xi);
  exprh = XISCALE*sinh(rh);


  xx = cosh(expxi)*cos(th)*cosh(exprh)*sin(ph) + sinh(expxi)*sin(th)*sinh(exprh)*cos(ph);
  yy = sinh(expxi)*sin(th)*cosh(exprh)*sin(ph) - cosh(expxi)*cos(th)*sinh(exprh)*cos(ph);

  rr = pow(xx*xx+yy*yy,1.0/Nval);
  angle = (2.0*atan2(yy,xx) +  M_PI + 2*K1val*M_PI)/Nval;

  X1 = rr*cos(angle);
  Y1 = rr*sin(angle);

  xx =  cosh(expxi)*sin(th)*cosh(exprh)*sin(ph) - sinh(expxi)*cos(th)*sinh(exprh)*cos(ph); 
  yy = -sinh(expxi)*cos(th)*cosh(exprh)*sin(ph) - cosh(expxi)*sin(th)*sinh(exprh)*cos(ph);

  rr = pow(xx*xx+yy*yy,1.0/Nval);
  angle = (2.0*atan2(yy,xx) +  M_PI + 2*K2val*M_PI)/Nval;

  X2 = rr*cos(angle);
  Y2 = rr*sin(angle);

  xx = cosh(exprh)*cos(ph);
  yy = sinh(exprh)*sin(ph);

  rr = pow(xx*xx+yy*yy,1.0/Nval);
  angle = (2.0*atan2(yy,xx) +  M_PI + 2*K3val*M_PI)/Nval;

  X3 = rr*cos(angle);
  Y3 = rr*sin(angle);


  X0=  1.0;  
  Y0 = 0.0;

  denom = (X0*X0 + Y0*Y0 + X1*X1 + Y1*Y1 + X2*X2 + Y2*Y2 +X3*X3 + Y3*Y3);
 
  DIM(0) =  (X0*X0 + Y0*Y0)/denom;
  DIM(1) =  (X1*X1 + Y1*Y1)/denom;
  DIM(2) =  (X2*X2 + Y2*Y2)/denom;
  DIM(3) =  (X3*X3 + Y3*Y3)/denom;
  DIM(4) =  rt2*(X0* X1 + Y0* Y1)/denom;
  DIM(5) =  rt2*(X1* Y0 - X0* Y1)/denom;
  DIM(6) =  rt2*(X0* X2 + Y0* Y2)/denom;
  DIM(7) =  rt2*(X0* Y2 - X2* Y0)/denom;
  DIM(8) =  rt2*(X0* X3 + Y0* Y3)/denom;
  DIM(9) =  rt2*(X0* Y3 - X3* Y0)/denom;
  DIM(10) =  rt2*(X2* X3 + Y2* Y3)/denom;
  DIM(11) = -  rt2*(X3* Y2 - X2* Y3)/denom; /* This minus sign appears necessary */
  DIM(12) =  rt2*(X3* X1 + Y3* Y1)/denom;
  DIM(13) =  rt2*(X1* Y3 - X3* Y1)/denom;
  DIM(14) =  rt2*(X1* X2 + Y1* Y2)/denom;
  DIM(15) =  rt2*(X2* Y1 - X1* Y2)/denom;

}




///////////////////////////////////////////////////////////////
// (PART C) DEFINE THE ND EQUATION
///////////////////////////////////////////////////////////////


// return false if this equation is not defined.

PLUGIN_API bool
defineEquation( ND_Equation_Def *eqdef )
{
    ///////////////////////////////////////////////////////
    // (1) Comment the first line below to when the code is ready
    ///////////////////////////////////////////////////////

    // *** <4> EDIT HERE

    //#define EQUATION_NOT_DEFINED

    #ifdef EQUATION_NOT_DEFINED

    return false;

    #else


    ///////////////////////////////////////////////////////
    // (2) Define the name and description
    ///////////////////////////////////////////////////////

    // *** <5> EDIT HERE

    eqdef->name        = strdup( "CP3 Kummer - 00"   ) ;
    eqdef->description = strdup( K3_000_DESC ) ;


    ///////////////////////////////////////////////////////
    // (3) Number of Parameters (M) and
    //     Number of Dimension  (N)
    ///////////////////////////////////////////////////////

    eqdef->numPara = NUM_PARA ;
    eqdef->numDim  = NUM_DIM  ;


    ///////////////////////////////////////////////////////
    // (4) Parameters Range and Max. Subdivision
    ///////////////////////////////////////////////////////

    // - Make sure the array index isn't out of range, otherwise, may crash!
    // - from 0 to NUM_PARA-1

    // *** <6> EDIT HERE    th = theta is first, then xi

    eqdef->range_min[0] = 0.0 ;
    eqdef->range_max[0] = M_PI/2  ; // theta 0 to pi/2  
    eqdef->max_subdivision[0] = THSTEPS ;

    eqdef->range_min[1] = 0.0;
      eqdef->range_max[1] = XIMAX ; // -ximax to ximax
      eqdef->max_subdivision[1] = XISTEPS ;

  //      eqdef->range_min[1] = 0.0 ;
  //      eqdef->range_max[1] = M_PI/2  ; // theta 0 to pi/2
  //     eqdef->max_subdivision[1] = THSTEPS ;

   //    eqdef->range_min[3] = - XIMAX;
   //    eqdef->range_max[3] = XIMAX ; // -ximax to ximax
   //    eqdef->max_subdivision[3] = XISTEPS ;




    return true;


    #endif
}




///////////////////////////////////////////////////////////////
// (PART D) FREE THE ALLOCATED MEMORY (DO NOT EDIT)
///////////////////////////////////////////////////////////////

PLUGIN_API void
freeEquation( const ND_Equation_Def *eqdef )
{
    if ( eqdef->description ) free ( eqdef->description ) ;
    if ( eqdef->name        ) free ( eqdef->name        ) ;
}

////////////////////////////////////////////////////
// (PART E) Optional function for setting the color function
////////////////////////////////////////////////////

// Note:
// - NUM_PARA - number of elements in in_para_array
// - there are three colors in the out_color_array (output: ranged [0,1])
// - Make sure the array index isn't out of range, otherwise, may crash!

// *** <7> MAY EDIT HERE

PLUGIN_API void
colorFunction( const double in_para_array [_MAX_PARAM] , double out_color_array [3] )
{
    out_color_array[0] = 0.5 ;
    out_color_array[1] = 0.5 ;
    out_color_array[2] = 1.0 ;
}