(*                                   ASSIGNMENT 3a 
                                    Mandelbrot set.
  How to use.
  -----------
     On the sml prompt type:
 
      - open "complex.mlc";

     To make the fractal: 
      Syntax of makefrac:
     fun makefrac(z,c1,c2 :COMPLEX,e,thresh :real,max :int) 
       where z, c1,c2 are self-explanatory, e is the size of the square  
       for which one value of c is taken, thresh is the threshold value of 
       |z|, and max is the max. no. of iterations.

	the output is stored in the file frac1.tst
Note : 0.044 is a good value for e
Note 2 : To keep the no. of characters=1 for a pt. the following 
substitutions have been made : 

10 = #"A"
11 = #"B"

.
.
.
 
and so on.
*) 	

use "filedemo.sml";
signature COMP=
sig
type COMPLEX
val re :COMPLEX ->real
val img :COMPLEX ->real
val modu : COMPLEX ->real
val arg : COMPLEX -> real
val make :real*real ->COMPLEX
val realtocomp : real ->COMPLEX
val polar : COMPLEX -> real*real
val sum :COMPLEX*COMPLEX ->COMPLEX
val conj :COMPLEX->COMPLEX
val diff :COMPLEX*COMPLEX ->COMPLEX
val equalto: COMPLEX*COMPLEX*real ->bool
val prod :COMPLEX*COMPLEX ->COMPLEX
val square: COMPLEX -> COMPLEX
val divide  :COMPLEX*COMPLEX ->COMPLEX
val reciprocal :COMPLEX ->COMPLEX
val sin : COMPLEX ->COMPLEX
val cos : COMPLEX ->COMPLEX
val tan :COMPLEX ->COMPLEX
val exp :COMPLEX ->COMPLEX
val log :COMPLEX ->COMPLEX
val pow :COMPLEX*COMPLEX ->COMPLEX
val sinh :COMPLEX->COMPLEX
val cosh :COMPLEX->COMPLEX
val tanh :COMPLEX -> COMPLEX
val asin :COMPLEX ->COMPLEX
val acos :COMPLEX ->COMPLEX
val atan : COMPLEX -> COMPLEX
val asinh :COMPLEX ->COMPLEX
val acosh :COMPLEX ->COMPLEX
val atanh :COMPLEX ->COMPLEX
end

structure COM : COMP =
struct
datatype COMPLEX = COMPL of real*real
fun re(COMPL(a,_))=a
fun img(COMPL(_,b))=b
fun modu(COMPL(a,b))=Math.sqrt(a*a+b*b)
fun arg(COMPL(a,b))=Math.atan(b/a)
fun make(a,b)=COMPL(a,b)
fun realtocomp(a)=COMPL(a,0.0)
fun polar(a)=(modu(a),arg(a))
fun sum(COMPL(a1,b1),COMPL(a2,b2))=COMPL(a1+a2,b1+b2)
fun conj(COMPL(a,b))=COMPL(a,~b)
fun diff(COMPL(a1,b1),COMPL(a2,b2))=COMPL(a1-a2,b1-b2)
fun equalto(a,b,e)=modu(diff(a,b))<e
fun prod(COMPL(a1,b1),COMPL(a2,b2))=COMPL(a1*a2-b1*b2,a2*b1+a1*b2)
fun square(a)=prod(a,a)
fun divide(a,b)=let
	val COMPL(x,y)=prod(a,conj(b))
	val pp = Math.pow(modu(b),2.0)
	in COMPL(x/pp,y/pp) end
fun reciprocal(a)=divide(COMPL(1.0,0.0),a)
fun sin(COMPL(a,b))=
COMPL(Math.sin(a)*Math.cosh(b),Math.cos(a)*Math.sinh(b))
fun cos(COMPL(a,b))=
COMPL(Math.cos(a)*Math.cosh(b),Math.sin(a)*Math.sinh(b))
fun tan(a)=divide(sin(a),cos(a))
fun exp(COMPL(a,b))=COMPL(Math.exp(a)*Math.cos(b),Math.exp(a)*Math.sin(b))
fun log(COMPL(a,b))=COMPL(Math.ln(modu(COMPL(a,b))),Math.atan(b/a))
fun pow(a,b)=exp(prod(b,log(a)))
fun sinh(COMPL(a,b))=
COMPL(Math.sinh(a)*Math.cos(b),Math.cosh(a)*Math.sin(b))
fun cosh(COMPL(a,b))=
COMPL(Math.cosh(a)*Math.cos(b),Math.sinh(a)*Math.sin(b))
fun tanh(a)=divide(sinh(a),cosh(a))
fun asin(a)=
prod(COMPL(0.0,~1.0),log(sum(prod(COMPL(0.0,1.0),a),pow(diff(COMPL(1.0,0.0),pow(a,COMPL(2.0,0.0))),COMPL(0.5,0.0)))))
fun acos(a)=
prod(COMPL(0.0,~1.0),log(sum(prod(COMPL(0.0,1.0),a),pow(diff(pow(a,COMPL(2.0,0.0)),COMPL(1.0,0.0)),COMPL(0.5,0.0)))))
fun atan(a)=
prod(COMPL(0.0,0.5),log(divide(sum(COMPL(0.0,1.0),a),diff(COMPL(0.0,1.0),a))))
fun asinh(a)=
log (sum(a,pow(sum(pow(a,COMPL(2.0,0.0)),COMPL(1.0,0.0)),COMPL(0.5,0.0))))
fun acosh(a)=
log (sum(a,pow(diff(pow(a,COMPL(2.0,0.0)),COMPL(1.0,0.0)),COMPL(0.5,0.0))))					
fun atanh(a)=
prod(COMPL(0.0,0.5),log(divide(sum(COMPL(1.0,0.0),a),diff(COMPL(1.0,0.0),a))))
end;

open COM
(*Iterate on z*)
fun newz(z :COMPLEX, c:COMPLEX)=sum(c,z)

(*Check for divergence*)
fun check(z :COMPLEX, thresh :real)=(modu(z)>thresh)

(*find for a single number.*)
fun go(z,c :COMPLEX,count:int,max :int,thresh:real)=if count>max then max 
else
if check(z,thresh) then count else
go(newz(z,c),c,count+1,max,thresh)

fun conv(x :int)=if x <10 then Int.toString(x) else Char.toString(chr(65+x-10))

(*String output *)
fun finditr(z,c 
:COMPLEX,max:int,thresh:real)=conv(go(z,c,0,max,thresh)) :string
(*Output for the whole row*)
fun doforrow(x0,x1,y,e,thresh :real,max :int, z :COM.COMPLEX)=
if x0>x1 then "\n" else 
finditr(z,make(x0,y),max,thresh)^doforrow(x0+e,x1,y,e,thresh,max,z)

(*Complete output*)
fun findit(c1,c2,z :COM.COMPLEX, e,thresh :real,max 
:int)=
if img(c1)>img(c2) then "\n" else 
doforrow(re(c1),re(c2),img(c1),e,thresh,max,z)^findit(make(re(c1),img(c1)+e),c2,z,e,thresh,max)
fun makefrac(z,c1,c2 :COMPLEX,e,thresh :real,max :int)=
writefile("frac1.tst",findit(c1,c2,z,e,thresh,max))

(* For the suggested values *)

fun Mandelbrot(z,c1,c2:COMPLEX)=makefrac(z,c1,c2,0.03,2.0,50)
fun mandelbrot()=Mandelbrot(make(0.0,0.0),make(~2.25,~1.5),make(00.75,1.5))