#include <cmath>

#include <cassert>

#include <stdlib.h>

#include "Complex.h"

#include "SectionIntegrals.h"

/*

* C version of r1mach and d1mach from core (netlib)

* specialized for IEEE arithmetic

*

* MACHINE CONSTANTS (s for single, d for double)

* {S|D}1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.

* {S|D}1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.

* {S|D}1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING.

* {S|D}1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING.

* {S|D}1MACH(5) = LOG10(B)

*/

// #define M_PI 3.14159265358979323846264338328

// #define M_LN2 0.693147180559945309417232121458

static double mach[5] = {2.2250738585072014e-308,

1.7976931348623157e+308,

1.1102230246251565e-16,

2.2204460492503131e-16,

3.0102999566398120e-01};

static const double s11r[21] = {

0.0448875760891932036595562553276,0.0278480909574822965157922173757,

0.00394490790249120295818107628687,-0.00157697939158619172562804651751,

-0.0000886578217796691901712579357311,0.0000301708056772263120428135787035,

9.521839632337438230089618156e-7,-3.00028307455805582080773625835e-7,

-6.14917009583473496433650831019e-9,1.85133588988085286010092653662e-9,

2.67848449041765751590373973224e-11,-7.82394575359355297437491915705e-12,

-8.44240072511090922609176843848e-14,2.41333276776166240844516922196e-14,

2.02015531985181413114834031833e-16,-5.68171271075270422851146478874e-17,

-3.80082421064644521052871349836e-19,1.05551739229841670238163200361e-19,

5.7758422925275435667221605993e-22,-1.58774695838716531303310462626e-22,

-7.24181766014636685673730787292e-25

};

static const double s11i[21] = {

0.100116671557942715638078149123,0.0429600096728215971268270800599,

-0.00799014859477407505770275088389,-0.000664114111384495427035329182866,

0.000240714510952202000864758517061,9.89085259369337382687437812294e-6,

-3.22040860178194578481012477174e-6,-8.08401148192350365282200249069e-8,

2.48351290049260966544658921605e-8,4.24154988067028660399867468349e-10,

-1.25611378629704490237955971836e-10,-1.56053077919196502557674988724e-12,

4.50565044006801278137904597946e-13,4.2641179237225098728291226479e-15,

-1.2084245714879456268965803807e-15,-9.01338537885038989528688031325e-18,

2.5180796700698002962991581923e-18,1.51955263898294940481729370636e-20,

-4.19737873024216866691628952458e-21,-2.092488792285595339755624521e-23,

5.72708467031136321701747126611e-24

};

static const double s12r[21] = {

-0.376145877558191778393359413441,0.0775244431850198578126067647425,

0.0120396593748540634695397747695,-0.00385683684390247509721340352427,

-0.000232359275790231209370627606991,0.0000697318379146209092637310696007,

2.32354473986257272021507575389e-6,-6.71692140309360615694979580992e-7,

-1.43946361256617673523038166877e-8,4.06087820907414336567714443732e-9,

6.10183339004616075548375321861e-11,-1.69196418769523832825063863136e-11,

-1.88669746820541798989965091628e-13,5.16473095452962111184823547686e-14,

4.45066881692009291504139737861e-16,-1.20625107617859803735741992452e-16,

-8.28193837331508300767103116139e-19,2.22680015825230528892642524445e-19,

1.24755889505424049389100515561e-21,-3.33254971913153176741833960484e-22,

-1.55307002839777371508497520751e-24

};

static const double s12i[21] = {

0.0527472790869782317601048210983,0.00823962722148093961886198320927,

-0.0205185842051817330153151013327,-0.00184683218270819613487368071941,

0.000569681886932212757533488372406,0.0000248774530818801164177266528608,

-7.31121019876580624171992432347e-6,-1.92744564223806538367454388776e-7,

5.49794278719049727550379096876e-8,9.78237385539447442446850072421e-10,

-2.7341624177723508216430132999e-10,-3.51839815887772323640101921381e-12,

9.68934411607055794052256859665e-13,9.45703963505047353201918875825e-15,

-2.57516976113400217760868402425e-15,-1.97419921753098238455550504742e-17,

5.32820017906655555903355375475e-18,3.29581793797656865402793252539e-20,

-8.83137325823594007269279476114e-21,-4.50279718100548728336329365981e-23,

1.19941679774924468309434420379e-23

};

static const double m12r[21] = {

0.148523151773238914750879360089,-0.0117856118001224048185631301904,

-0.00248887208039014371691400683052,0.000250045060357076469386198883676,

0.0000227217776065076434637230864113,-2.48764935230787745662127026799e-6,

-1.32138506847814502856384193414e-7,1.50966754393693942843767293542e-8,

5.3472999553162661403204445045e-10,-6.26136041009708550772228055719e-11,

-1.59574066624737000616598104732e-12,1.89788785691219687197167013023e-13,

3.66030609080549274006207730375e-15,-4.39955659500182569051978906011e-16,

-6.65848768159000092224193226014e-18,8.06343127453005031535923212263e-19,

9.84397490339224661524630997726e-21,-1.19869887155210161836484730378e-21,

-1.20634550494837590549640883469e-23,1.47512193662595435067359954287e-24,

1.24549093756962710863096766634e-26

};

static const double m12i[21] = {

-0.0454399665519585306943416687117,-0.0210517666740874019203591488894,

0.00194647501081621201871675259482,0.000253466068123907163346571754613,

-0.0000268083453427538717591876419304,-1.82138740336918117478832696004e-6,

2.04357511048425337951376869602e-7,8.75944656915074206478854298947e-9,

-1.01466837126303146739791005703e-9,-3.02573132377805421636557302451e-11,

3.57358222114420372764650037191e-12,7.88121312149152771558608913996e-14,

-9.42758576193708862552405242331e-15,-1.60439904050827900099939709069e-16,

1.93624791035947590366500765061e-17,2.62394448214143482490534256935e-19,

-3.18700789496399461681365308408e-20,-3.52400207248027768109209530864e-22,

4.30074555255053206057921088056e-23,3.95655079023456015736315286131e-25,

-4.84642137915095135859812028886e-26

};

/*

* evaluate a chebyshev series

* adapted from fortran csevl

*/

static double

csevl( const double x, const double *cs, unsigned int n )

{

double b2 = 0, b1 = 0, b0 = 0;

const double twox = 2 * x;

while( n-- ){

b2 = b1;

b1 = b0;

b0 = twox * b1 - b2 + cs[n];

}

return ( b0 - b2 )/2;

}

/*

* from the original fortran inits

* april 1977 version. w. fullerton, c3, los alamos scientific lab.

*

* initialize the orthogonal series so that inits is the number of terms

* needed to insure the error is no larger than eta. ordinarily, eta

* will be chosen to be one-tenth machine precision.

*/

static unsigned int

inits( const double *series, unsigned int n, const double eta ){

double err = 0;

while( err <= eta && n-- ){

err += fabs( series[n] );

}

return n++;

}

static DDG::Complex

s11( const double t ){

// how many entries to use?

static unsigned int ns11r = 0, ns11i = 0;

if( !ns11r ){

ns11r = inits( s11r, sizeof s11r / sizeof *s11r, mach[2] / 10 );

ns11i = inits( s11i, sizeof s11i / sizeof *s11i, mach[2] / 10 );

}

return DDG::Complex( csevl( t, s11r, ns11r ), csevl( t, s11i, ns11i ) );

}

static DDG::Complex

s12( const double t ){

// how many entries to use?

static unsigned int ns12r = 0, ns12i = 0;

if( !ns12r ){

ns12r = inits( s12r, sizeof s12r / sizeof *s12r, mach[2] / 10 );

ns12i = inits( s12i, sizeof s12i / sizeof *s12i, mach[2] / 10 );

}

return DDG::Complex( csevl( t, s12r, ns12r ), csevl( t, s12i, ns12i ) );

}

static DDG::Complex

m12( const double t ){

// how many entries to use?

static unsigned int nm12r = 0, nm12i = 0;

if( !nm12r ){

nm12r = inits( m12r, sizeof m12r / sizeof *m12r, mach[2] / 10 );

nm12i = inits( m12i, sizeof m12i / sizeof *m12i, mach[2] / 10 );

}

return DDG::Complex( csevl( t, m12r, nm12r ), csevl( t, m12i, nm12i ) );

}

DDG::Complex

DirichletIJDirect( const double s, const double gii, const double gij, const double gjj ){

const double s2(s*s);

const DDG::Complex is(0,s), is3(0,s*s2);

const double s4(s2*s2);

return

(((3*gii + 4*gij + 3*gjj) + is*(gii + gij + gjj) - is3*gij/6

+ DDG::Complex(cos(s),sin(s))*(-(3*gii + 4*gij + 3*gjj)

+ is*(2*gii + 3*gij + 2*gjj)

+ s2*(gii + 2*gij + gjj)/2))/s4

+ (gii - 2*gij + gjj)/24 - is*(gii - 2*gij + gjj)/60);

}

// gii = |e_{ki}|^2, gij = <e_{ki}, e_{kj}>, gjj = |e_{jk}|^2

DDG::Complex

DirichletIJ( const double s, const double gii, const double gij, const double gjj ){

if( fabs(s) > M_PI ){

return DirichletIJDirect( s, gii, gij, gjj );

}else if( s > 0 ){

// between 0 and pi

const double t = s*2/M_PI - 1; // normalize

return ((gii + gjj)*s11(t) + gij*s12(t)).conj();

}else{

// between -pi and 0

const double t = -s*2/M_PI - 1; // normalize

return ((gii + gjj)*s11(t) + gij*s12(t));

}

}

// gjj = |e_{ij}|^2, gjk = <e_{ij}, e_{ik}>, gkk = |e_{ki}|^2

double

DirichletII( const double s, const double gjj, const double gjk, const double gkk ){

return .25*((gjj-2*gjk+gkk)+s*s*(gjj+gjk+gkk)/90);

}

double

MassII( void ){ return 1/6.; }

DDG::Complex

MassIJDirect( const double s ){

const double s2 = s*s;

const double s4 = s2*s2;

const DDG::Complex is(0,s), is3(0,s*s2);

return (6*DDG::Complex(cos(s), sin(s)) - 6 - 6*is + 3*s2 + is3)/(3*s4);

}

DDG::Complex

MassIJ( const double s ){

if( fabs(s) > M_PI ){

return MassIJDirect( s );

}else if( s > 0 ){

const double t = s/M_PI*2-1;

return m12(t).conj();

}else{

const double t = -s/M_PI*2-1;

return m12(t);

}

}