// This program illustrates using user-defined functions to perform
// complex math and trig functions, logarithms, etc.
// Since user-defined functions can only return one value, 
//   separate functions return the real and imaginary parts
//   of each function.
// Refer to a textbook on Complex Variables.

long "all"	// not required, unless long precision is desired.

100
input "Enter real, imaginary: "; in_real, in_imag
// Keep in mind the range of each function when entering values.
// (small numbers work best)

// calculate exponential
exp_real = fnexpz_real (in_real, in_imag)
exp_imag = fnexpz_imag (in_real, in_imag)
print "exp(z) [real]:   "; exp_real
print "       [imag]:   "; exp_imag

// calculate log
log_real = fnlogz_real (exp_real, exp_imag)
log_imag = fnlogz_imag (exp_real, exp_imag)
print "log(z) [real]:   "; log_real
print "       [imag]:   "; log_imag

// calculate sin
sin_real = fnsinz_real (in_real, in_imag)
sin_imag = fnsinz_imag (in_real, in_imag)
print "sin(z) [real]:   "; sin_real
print "       [imag]:   "; sin_imag

// calculate cosine
cos_real = fncosz_real (in_real, in_imag)
cos_imag = fncosz_imag (in_real, in_imag)
print "cos(z) [real]:   "; cos_real
print "       [imag]:   "; cos_imag

// Calculate sin^2 + cos^2  (should equal 1)
sx2 = fnmultz_real (sin_real, sin_imag, sin_real, sin_imag)
sy2 = fnmultz_imag (sin_real, sin_imag, sin_real, sin_imag)
cx2 = fnmultz_real (cos_real, cos_imag, cos_real, cos_imag)
cy2 = fnmultz_imag (cos_real, cos_imag, cos_real, cos_imag)
real = sx2 + cx2
imag = sy2 + cy2
print "s2+c2  [real]:   "; real
print "       [imag]:   "; imag

// Calculate tangent
tan_real = fntanz_real (in_real, in_imag)
tan_imag = fntanz_imag (in_real, in_imag)
print "tan    [real]:   "; tan_real
print "       [imag]:   "; tan_imag
 
// Calculate arctangent - should yield original values
//                        if range is not exceeded.
atn_real = fnatnz_real (tan_real, tan_imag)
atn_imag = fnatnz_imag (tan_real, tan_imag)
print "atn    [real]:   "; atn_real
print "       [imag]:   "; atn_imag

// Calculate cube using log and exp
if (in_real = 0) then goto 300  // win't work for pure imaginary
log3_real = fnlogz_real (in_real, in_imag)
log3_imag = fnlogz_imag (in_real, in_imag)
log3_3_real = fnmultz_real (log3_real, log3_imag, 3, 0)
log3_3_imag = fnmultz_imag (log3_real, log3_imag, 3, 0)
cube_real = fnexpz_real (log3_3_real, log3_3_imag)
cube_imag = fnexpz_imag (log3_3_real, log3_3_imag)
print "cube   [real]:   "; cube_real
print "       [imag]:   "; cube_imag

// Now calculate the cube root - should return original number
log3_real = fnlogz_real (cube_real, cube_imag)
log3_imag = fnlogz_imag (cube_real, cube_imag)
log3_3_real = fndivz_real (log3_real, log3_imag, 3, 0)
log3_3_imag = fndivz_imag (log3_real, log3_imag, 3, 0)
cuberoot_real = fnexpz_real (log3_3_real, log3_3_imag)
cuberoot_imag = fnexpz_imag (log3_3_real, log3_3_imag)
print "cube-  [real]:   "; cuberoot_real
print "root   [imag]:   "; cuberoot_imag
300

// Now calculate the exp (i * pi) -- should equal -1
eipi_real = fnexpz_real (0, pi)
eipi_imag = fnexpz_imag (0, pi)
print "e^ipi  [real]:   "; eipi_real
print "       [imag]:   "; eipi_imag

// Calculate the square root
sqrt_real = fnsqrz_real (in_real, in_imag)
sqrt_imag = fnsqrz_imag (in_real, in_imag)
print "sqrt   [real]:   "; sqrt_real
print "       [imag]:   "; sqrt_imag

goto 100

//  USER-DEFINED FUNCTIONS FOR COMPLEX OPERATIONS:
//--------------------------------------
// Complex Exponential -- exp (x + iy)
def fnexpz_real(x, y) = exp(x) * cos (y)
def fnexpz_imag(x, y) = exp(x) * sin (y)
//-----------------------------------------
// Complex natural log -- ln (x + iy)
// works for primary complex plane (-pi < theta < pi)
def fnlogz_real (x,y) = log (sqr(x*x + y*y))
def fnlogz_imag (x,y) = gosub 1250
1250 // complex log
// works for primary complex plane (-pi < theta < pi)
// if x = 0, we have a problem
if (x < 0 and y >= 0) then return atn(y/x) + pi
if (x < 0 and y < 0)  then return atn(y/x) - pi
return atn (y/x)
//---------------------------------------------
// Hyperbolic sine and cosine (needed for complex sine and cosine)
def fnsinh (y) = (exp(y) - exp(-y))/2
def fncosh (y) = (exp(y) + exp(-y))/2
//---------------------------------------------
// Complex sine -- sin (x + iy)
def fnsinz_real (x, y) = sin (x) * fncosh (y)
def fnsinz_imag (x, y) = cos (x) * fnsinh (y)
//---------------------------------------------
// Complex cosine -- cos (x + iy)
def fncosz_real (x, y) = cos  (x) * fncosh (y)
def fncosz_imag (x, y) = -sin (x) * fnsinh (y)
//---------------------------------------------
// complex multiply -- (x1 + iy1) * (x2 + iy2)
def fnmultz_real (x1, y1, x2, y2)= x1 * x2 - y1 * y2
def fnmultz_imag (x1, y1, x2, y2)= x1 * y2 + x2 * y1
//---------------------------------------------
// complex divide -- (x1 + iy1) / (x2 + iy2)
def fndivz_real  (x1, y1, x2, y2) = (x1 * x2 + y1 * y2) / (x2 * x2 + y2 * y2)
def fndivz_imag  (x1, y1, x2, y2) = (y1 * x2 - x1 * y2) / (x2 * x2 + y2 * y2)
//----------------------------------------
// Complex tangent -- tan (x + iy)
def fntanz_real (x, y) = fndivz_real (fnsinz_real(x,y), fnsinz_imag(x,y), fncosz_real (x,y), fncosz_imag (x,y))
def fntanz_imag (x, y) = fndivz_imag (fnsinz_real(x,y), fnsinz_imag(x,y), fncosz_real (x,y), fncosz_imag (x,y))
//----------------------------------------
// Complex arctangent -- atn (x + iy)
// atn (z) = 1/2 i (ln(1-iz) - ln(1+iz))
//	   = 1/2 i (ln (1 - i (x + iy))   -   ln (1 + i (x + iy))
//         = 1/2 i (ln (1 + y, -ix)       -   ln (1 - y + ix))
def fnatnz_real (x, y) = 0.5 * (fnlogz_imag (1-y, x) - fnlogz_imag (1+y, -x))
def fnatnz_imag (x, y) = 0.5 * (fnlogz_real (1+y, -x) - fnlogz_real (1-y, x))
//-----------------------------------------------
// Complex square root -- sqr (x + iy)
def fnsqrz_real(x,y) = sqr ((x + sqr (x*x + y*y))/ 2)
def fnsqrz_imag(x,y) = (sgn(y)+(y = 0)) * sqr((sqr(x*x + y*y) - x)/2)
// Note: (sgn(y) + (y = 0)) gives sgn(y) or 1 if y=0
//-----------------------------------------------