#################################################
#### COLLECTION OF TOOLS TO WORK WITH MOEBIUS TRANSFORMATIONS #######
#### © MIKA SEPPÄLÄ 2002
#################################################


CR := proc(a,b,c,d)
        if b=c or c=d or b=d then
            ERROR("Cross ratio not defined.")
        fi;
        if a<>  infinity then
        if d = infinity then
           (a-c)/(b-c);
        elif c = infinity then
           (b-d)/(a-d);
        elif b = infinity then
           (a-c)/(a-d);
        else
           (a-c)*(b-d)/(a-d)/(b-c);
        fi;
        else
        if d = infinity then
            infinity;
        elif c = infinity then
            0;
        elif b = infinity then
           1;
        else
          (b-d)/(b-c);
        fi;
        fi;
  end;

Mat2Moe := proc(M, matrix)
       if linalg[det](M)=0 then ERROR("Give a matrix with non-zero determinant."); fi;
       unapply((M[1,1]*z + M[1,2])/(M[2,1]*z + M[2,2]),z);
  end;

Moe2Mat := proc(g)
      local M, nu,de,i,j;
      nu := collect(numer(simplify(g(z))),z): de := collect(denom(simplify(g(z))),z):
      if degree(nu,z)>1 then ERROR("Input not a Moebius transformation"); fi;
      if degree(de,z)>1 then ERROR("Input not a Moebius transformation"); fi;
      M := matrix(2,2);
      for i from 1 to 2 do
          M[1,i] := coeff(nu,z,2-i):
      od;
      for j from 1 to 2 do
          M[2,j] := coeff(de,z,2-j):
      od;
         if linalg[det](evalm(M)) = 0 then
           ERROR("Give a non-constant mapping");
         fi;
  if M[1,1] + M[2,2]  = 0 then map(x -  x/sqrt(linalg[det](M)), evalm(M));
  else
  map(x-> -x/sqrt(linalg[det](M)), evalm(M)); fi;
  end;

 CreateEllMoe := proc(a,r,theta)
  if a=r then
     ERROR("Fix points (the first two argument) must be distinct.");     fi;
  if theta=infinity then
     ERROR("The last argument (rotation angle) must be between 0 and 2*Pi ");
  fi;
#  if not(theta =0) or not(theta =evalf(2*Pi)) then
#     ERROR("The last argument (rotation angle) must be between 0 and #2*Pi.");
#  fi;
  if a <  infinity and r <  infinity then
    unapply(((exp(I*theta)*a-r)*z-a*r*(exp(I*theta)-1))/((exp(I*theta)-1)*z + a - exp(I*theta)*r),z);
  elif a = infinity then
    unapply(exp(I*theta)*z - r*(exp(I*theta)-1),z);
  else
    unapply(z/exp(I*theta) + a*(1-1/exp(I*theta)),z);
  fi;
end;
CreateHypMoe := proc(a,r,k)
    if a=r then
       ERROR("Fix points (the first two argument) must be distinct.");     fi;
    if k=infinity then
       ERROR("The last argument (multiplier) must be a positive number.");
    fi;
    if not(k > 0) then
       ERROR("The last argument (multiplier) must be a positive number.");
    fi;
    if a <  infinity and r <  infinity then
      unapply(((k*a-r)*z-a*r*(k-1))/((k-1)*z + a - k*r),z);
    elif a = infinity then
      unapply(k*z - r*(k-1),z);
    else
      unapply(z/k + a*(1-1/k),z);
    fi;
  end;


  HyperMoeFixPoints := proc(g)
        local x,y,a,r;
        if WhatType(g) <  hyperbolic then
           ERROR("Not a hyperbolic Moebius transformation.");
        fi;
        if abs(limit(g(z),z=infinity))=infinity then
           x:= infinity;
           y := solve(g(z)=z,z);
         if abs(g(y+I)-y)<1 then
           a,r := y,x;
         else
           a,r := x,y;
         fi;
       else
       x := solve(g(z)=z,z)[1]; y := solve(g(z)=z,z)[2];
       if abs(g((x+y)/2)-x) < abs((y-x)/2) then
           a,r := x,y;
       else
           a,r := y,x;
       fi;
       fi;
       a,r;
  end;

  Multiplier := proc(g)
       local a,r;
       a,r := HyperMoeFixPoints(g);
       if a = infinity then
       CR(g(r+1),r+1,r,a);
       elif r = infinity then
       CR(g(a+1),a+1,r,a);
       else
       CR(g((a+r)/2),(a+r)/2,r,a);
       fi;
  end;



WhatType := proc(g)
        local M;
        if diff(g(z),z)=1 then
            identity;
        else
        M := Moe2Mat(g);
        if evalf(abs(linalg[trace](M)))> 2 then
            hyperbolic;
        elif evalf(abs(linalg[trace](M)))=2 then
            parabolic;
        else
            elliptic;
        fi;
        fi;
  end;

WhatTypeM := proc(g)
       local x,y, TYPE;
       if diff(g(z),z)=1 then
           identity;
       else
       if Im(g(I))<=0 then
         ERROR("Tansformation does not map the upper half-plane onto itself.");
       fi;
       if abs(limit(g(z),z=infinity))=infinity then
           if nops([solve(g(z)=z,z)])=0 then TYPE := 'parabolic';
           elif nops([solve(g(z)=z,z)])=1 and type(solve(g(z)=z,z),numeric) then
             TYPE := 'hyperbolic';
           elif nops([solve(g(z)=z,z)])=1 and not(type(solve(g(z)=z,z),numeric)) then
             TYPE := 'parabolic'
           else ERROR("Wrong type of input.")
       fi;
       fi;
       if abs(limit(g(z),z=infinity))< infinity
       then if
           nops([solve(g(z)=z,z)])=2 and Im(solve(g(z)=z,z)[1])=0 then TYPE := 'hyperbolic';
       elif
           nops([solve(g(z)=z,z)])=2 and Im(solve(g(z)=z,z)[1])< 0 then TYPE := 'elliptic';
       elif
             nops([solve(g(z)=z,z)])=1 then TYPE := 'parabolic'; fi;
       fi;
       fi;
       TYPE;
  end;