// complex_math.bal     type complex and many functions
//                      cx, real, imag, cxabs, cxconj, cxneg, arg, cxti,
//                      cxprtln, cxprt, cxsprtln,
//                      cxexp, cxlog, cssqrt,
//                      cxadd, cxsub, cxmul, cxdiv, 
//                      cxsin, cxcos, cxtan, cxasin, cxacos, csatan
//                      cxsinh, cxcosh, cxtanh,

import ballerina/io;
import ballerina/math;

type complex record {  // type complex  z = re + i*im
  float re;
  float im;
};

function cx (float x, float y) returns (complex) { // ceate z = x + iy
  complex z;
  z = {re: x, im: y}; // z = x + iy
  return z;
}

function real (complex z) returns (float) { // return real part
  return z.re;
}

function imag (complex z) returns (float) { // return imaginary part
  return z.im;
}

function cxabs (complex z) returns (float) { // absolute value, also length
  return math:sqrt(z.re*z.re + z.im*z.im);
}

function cxconj (complex z) returns (complex) { // conjugate
  return cx(z.re, -z.im);
}

function cxneg (complex z) returns (complex) { // negate
  return cx(-z.re, -z.im);
}

function cxarg (complex z) returns (float) { // argument, angle
  return math:atan2(z.im, z.re);
}

function cxti (complex z) returns (complex) { // times i = i * z
  return cx(-z.im, z.re);
}

function cxprtln (complex z) { // print (x,y) format
  float x = z.re;
  float y = z.im;
  io:println("("+x+","+y+")");
}

function cxprt (complex z) { // print, no line feed  (x,y) format
  float x = z.re;
  float y = z.im;
  io:print("("+x+","+y+")");
}

function cxsprtln(string s, complex z) { // print s = z
  float x = z.re;
  float y = z.im;
  io:println(s+" = ("+x+","+y+")");
}

function cxadd (complex z1, complex z2) returns (complex) { // add 
  return {re: (z1.re+z2.re), im: (z1.im+z2.im)};
}

function cxsub (complex z1, complex z2) returns (complex) { // subtract
  return {re: (z1.re-z2.re), im: (z1.im-z2.im)};
}

function cxmul (complex z1, complex z2) returns (complex) { // multiply
  return {re: (z1.re*z2.re - z1.im*z2.im),
          im: (z1.re*z2.im + z1.im*z2.re)};
}

function cxdiv (complex z1, complex z2) returns (complex) { // divide  
  float r = z2.re*z2.re + z2.im*z2.im;
  return {re: (z1.re*z2.re + z1.im*z2.im)/r,
          im: (z1.im*z2.re - z1.re*z2.im)/r};
}

function cxexp (complex z) returns (complex) { // exp
  float expx = math:exp(z.re);
  return cx(expx*math:cos(z.im),
            expx*math:sin(z.im));
}

function cxlog (complex z) returns (complex) { // log
  float PI = 3.141592653589793238462643;
  float rpart = math:log(math:sqrt(z.re*z.re+z.im*z.im));
  float ipart = math:atan2(z.im, z.re);
  if (ipart>PI) { ipart = ipart-2.0*PI; }
  return cx(rpart, ipart);
}

function cxsqrt (complex z) returns (complex) { // square root
  float r = math:sqrt(z.re*z.re + z.im*z.im);
  float rpart = math:sqrt(0.5*(r+z.re));
  float ipart;
  ipart = math:sqrt(0.5*(r-z.re));
  if (z.im<0.0) { ipart = -ipart; }
  return cx(rpart, ipart);
}

function cxsin (complex z) returns (complex) { // sine
  return cx(math:sin(z.re)*math:cosh(z.im),
            math:cos(z.re)*math:sinh(z.im));
}

function cxcos (complex z) returns (complex) { // cosine
  return cx(math:cos(z.re)*math:cosh(z.im),
            -math:sin(z.re)*math:sinh(z.im));
}

function cxtan (complex z) returns (complex) { // tangent
  return cxdiv(cxsin(z),cxcos(z));
}

function cxcot (complex z) returns (complex) { // co tangent
  return cxdiv(cxcos(z),cxsin(z));
}

function cxasin (complex z) returns (complex) { // arcsine
  // -i*log(i*z + sqrt(1-z)*sqrt(1+z) )
  complex z1 = cxsqrt(cxsub(cx(1.0,0.0),z));
  complex z2 = cxsqrt(cxadd(cx(1.0,0.0),z));
  complex z3 = cxadd(cxti(z),cxmul(z1,z2));
  return cxneg(cxti(cxlog(z3)));
}

function cxacos (complex z) returns (complex) { // arccosine
  // -i * log(z + sqrt(1-z^2) )
  complex zz = cxmul(z,z);
  complex sq = cxsqrt(cxsub(cx(1.0,0.0), zz));
  complex lz = cxlog(cxadd(z, sq));
  return cxneg(cxti(lz));
}

function cxatan (complex z) returns (complex) { // arctangent
  complex z1 = cx(0.0, -1.0);
  complex z2 = cx(z.re, z.im-1.0);
  complex z3 = cx(-z.re, -z.im-1.0);
  return cx(0.0, 0.0);
}

public function main (string... args) {
  io:println("complex_math.bal running");

  complex z1;
  complex z2;
  complex z3;
  complex z4;
  complex zt;
  float x;
  float y;
  float a;
  float err;
  
  z1 = {re: 2.0, im: 3.0};  
  io:print("z1 = ");
  cxprtln(z1);
  z2 = cx(-2.0, 3.0);
  cxprt(z2);
  io:println(" = z2");
  cxsprtln("z2",z2);
  z3 = cx(2.0, -3.0);
  cxsprtln("z3",z3);
  z4 = cx(-2.0, -3.0); // four test values with different signs
  cxsprtln("z4",z4);

  io:println(" ");
  x = real(z1);
  io:println("x=real(z1) = "+x);
  y = imag(z1);
  io:println("y=imag(z1) = "+y);
  a = cxabs(z1);
  io:println("a=cxabs(z1) = "+a);
  zt = cxconj(z1);
  cxsprtln("zt=cxconj(z1)",zt);
  zt = cxneg(z1);
  cxsprtln("zt=cxneg(z1)",zt);
  a = cxarg(z1);
  io:println("a=cxarg(z1) = "+a);
  zt = cxti(z1);
  cxsprtln("zt=cxti(z1)",zt);

  io:println(" ");
  zt = cxadd(z1, z2);
  cxsprtln("zt=z+z2",zt);
  zt = cxsub(z1, z2);
  cxsprtln("zt=z-z2",zt);
  zt = cxmul(z1, z2);
  cxsprtln("zt=z*z2",zt);
  zt = cxdiv(z1, z2);
  cxsprtln("zt=z/z2",zt);
  z4 = cxmul(zt, z2);
  cxsprtln("z4=zt*z2=z? ",z4);

  io:println(" ");
  zt = cxexp(z1);
  cxsprtln("zt = cxexp(z1)",zt);
  zt = cxlog(z1);
  cxsprtln("zt = cxlog(z1)",zt);
  zt = cxsqrt(z1);
  cxsprtln("zt = cxqrt(z1)",zt);

  io:println(" ");
  zt = cxsin(z1);
  cxsprtln("zt = cxsin(z1)",zt);
  zt = cxcos(z1);
  cxsprtln("zt = cxcos(z1)",zt);
  zt = cxtan(z1);
  cxsprtln("zt = cxtan(z1)",zt);
  zt = cxcot(z1);
  cxsprtln("zt = cxcot(z1)",zt);
  io:println(" ");
  zt = cxasin(z1);
  cxsprtln("zt = cxasin(z1)",zt);
  zt = cxacos(z1);
  cxsprtln("zt = cxacos(z1)",zt);
  zt = cxasin(z1);
  cxsprtln("zt = cxatan(z1)",zt);
  io:println(" ");

  io:println("test to check accuracy");
  zt = cxsub(cxadd(z1,z2),z2); 
  err = cxabs(cxsub(z1,zt));
  io:println("err zt-cxsub(cxadd(z1,z2),z2) = "+err);
 
  zt = cxdiv(cxmul(z1,z2),z2); 
  err = cxabs(cxsub(z1,zt));
  io:println("err zt-cxdiv(cxmul(z1,z2),z2) = "+err);
 
  zt = cxmul(cxsqrt(z1),cxsqrt(z1)); 
  cxsprtln("cxsqrt(z1)*cxsqrt(z1)=z1",zt);
  err = cxabs(cxsub(z1,zt));
  io:println("err zt-cxsqrt(z1)*cxsqrt(z1) = "+err);
  zt = cxmul(cxsqrt(z2),cxsqrt(z2)); 
  cxsprtln("cxsqrt(z2)*cxsqrt(z2)=z2",zt);
  err = cxabs(cxsub(z2,zt));
  io:println("err zt-cxsqrt(z1)*cxsqrt(z2) = "+err);
  zt = cxmul(cxsqrt(z3),cxsqrt(z3)); 
  cxsprtln("cxsqrt(z3)*cxsqrt(z3)=z3",zt);
  err = cxabs(cxsub(z3,zt));
  io:println("err zt-cxsqrt(z3)*cxsqrt(z3) = "+err);
  zt = cxmul(cxsqrt(z4),cxsqrt(z4)); 
  cxsprtln("cxsqrt(z4)*cxsqrt(z4)=z4",zt);
  err = cxabs(cxsub(z4,zt));
  io:println("err zt-cxsqrt(z4)*cxsqrt(z4) = "+err);
 
  zt = cxlog(cxexp(z1)); 
  cxsprtln("cxlog(cxexp(z1))=z1",zt);
  err = cxabs(cxsub(z1,zt));
  io:println("err zt-cxlog(cxexp(z1)) = "+err);
 
  zt = cxlog(cxexp(z2)); 
  cxsprtln("cxlog(cxexp(z2))=z2",zt);
  err = cxabs(cxsub(z2,zt));
  io:println("err zt-cxlog(cxexp(z2)) = "+err);
 
  zt = cxexp(cxlog(z1)); 
  cxsprtln("cxexp(cxlog(z1))=z1",zt);
  err = cxabs(cxsub(z1,zt));
  io:println("err zt-cxexp(cxlog(z1)) = "+err);
 
  zt = cxexp(cxlog(z2)); 
  err = cxabs(cxsub(z2,zt));
  io:println("err z2-cxexp(cxlog(z2)) = "+err);

  zt = cxdiv(cxsub(cxexp(cxti(z1)),cxexp(cxneg(cxti(z1)))),cx(0.0,2.0));
  cxsprtln("(exp(i*z1)-exp(-i*z1)/2i=sin(z1)",zt);
  err = cxabs(cxsub(cxsin(z1),zt));
  io:println("err (exp(i*z1)-exp(-i*z1))/2i-sin(z1)) = "+err);

  zt = cxdiv(cxadd(cxexp(cxti(z1)),cxexp(cxneg(cxti(z1)))),cx(2.0,0.0));
  cxsprtln("(exp(i*z1)+exp(-i*z1)/2=cos(z1)",zt);
  err = cxabs(cxsub(cxcos(z1),zt));
  io:println("err (exp(i*z1)+exp(-i*z1))/2-cos(z1)) = "+err);

  zt = cxasin(cxsin(z1));
  cxsprtln("cxasin(cxsin(z1))=z1",zt);
  err = cxabs(cxsub(z1,zt));
  io:println("err cxasin(cxsin(z1)) = "+err);
  zt = cxsin(cxasin(z1));
  err = cxabs(cxsub(z1,zt));
  io:println("err cxsin(cxasin(z1)) = "+err);

  zt = cxacos(cxcos(z1));
  cxsprtln("cxacos(cxcos(z1))=z1",zt);
  err = cxabs(cxsub(z1,zt));
  io:println("err cxacos(cxcos(z1)) = "+err);
  zt = cxcos(cxacos(z1));
  err = cxabs(cxsub(z1,zt));
  io:println("err cxacos(cxacos(z1)) = "+err);

  zt = cxadd(cxmul(cxsin(z1),cxsin(z1)),
             cxmul(cxcos(z1),cxcos(z1)));
  cxsprtln("on z1 sin^2+cos^2",zt);
  zt = cxadd(cxmul(cxsin(z2),cxsin(z2)),
             cxmul(cxcos(z2),cxcos(z2)));
  cxsprtln("on z2 sin^2+cos^2",zt);
  zt = cxadd(cxmul(cxsin(z3),cxsin(z3)),
             cxmul(cxcos(z3),cxcos(z3)));
  cxsprtln("on z3 sin^2+cos^2",zt);
  zt = cxadd(cxmul(cxsin(z4),cxsin(z4)),
             cxmul(cxcos(z4),cxcos(z4)));
  cxsprtln("on z4 sin^2+cos^2",zt);
	     

  io:println(" ");
  io:println("complex_math.bal finished");
} // end complex_math.bal