#include "functions.inc"

#macro DefineOrbital(n,l,m,Z,value)

  #local Fact = function(num) { prod(i,1,num,i) }
  #local a = abs(m);
  #local b = Z/n;
  #if (m = 0)
    #local n_Ang = sqrt((2*l+1)*Fact(l-a)/(4*pi*Fact(l+a)));
  #else
    #local n_Ang = sqrt((2*l+1)*Fact(l-a)/(2*pi*Fact(l+a)));
  #end
  #local n_Rad = sqrt(Fact(n-l-1)/(2*n*pow(Fact(n+l),3)))*pow(2*b,l+3/2)*Fact(n+l);

  #switch (l)
    #case (0)
      #declare f_Ang = function(x,y,z) { 1 }
    #break
    #case (1)
      #switch (m)
        #case (0)
          #declare f_Ang = function(x,y,z) { y }
        #break
        #case (1)
          #declare f_Ang = function(x,y,z) { z }
        #break
        #case (-1)
          #declare f_Ang = function(x,y,z) { x }
        #break
      #end
    #break
    #case (2)
      #switch (m)
        #case (0)
          #local n_Ang = n_Ang/2;
          #declare f_Ang = function(x,y,z) { 2*y*y-x*x-z*z }
        #break
        #case (1)
          #local n_Ang = n_Ang*3;
          #declare f_Ang = function(x,y,z) { y*z }
        #break
        #case (-1)
          #local n_Ang = n_Ang*3;
          #declare f_Ang = function(x,y,z) { y*x }
        #break
        #case (2)
          #local n_Ang = n_Ang*3;
          #declare f_Ang = function(x,y,z) { z*z-x*x }
        #break
        #case (-2)
          #local n_Ang = n_Ang*6;
          #declare f_Ang = function(x,y,z) { x*z }
        #break
      #end
    #break
    #case (3)
      #switch (m)
        #case (0)
          #local n_Ang = n_Ang/2;
          #declare f_Ang = function(x,y,z) { y*(2*y*y-3*x*x-3*z*z) }
        #break
        #case (1)
          #local n_Ang = n_Ang*3/2;
          #declare f_Ang = function(x,y,z) { (4*y*y-x*x-z*z)*z }
        #break
        #case (-1)
          #local n_Ang = n_Ang*3/2;
          #declare f_Ang = function(x,y,z) { (4*y*y-x*x-z*z)*x }
        #break
        #case (2)
          #local n_Ang = n_Ang*15;
          #declare f_Ang = function(x,y,z) { y*(z*z-x*x) }
        #break
        #case (-2)
          #local n_Ang = n_Ang*30;
          #declare f_Ang = function(x,y,z) { y*x*z }
        #break
        #case (3)
          #local n_Ang = n_Ang*15;
          #declare f_Ang = function(x,y,z) { z*(z*z-3*x*x) }
        #break
        #case (-3)
          #local n_Ang = n_Ang*15;
          #declare f_Ang = function(x,y,z) { x*(3*z*z-x*x) }
        #break
      #end
    #break
  #end

  #declare rho = function(x,y,z) { b*f_r(x,y,z) }
  #switch (n)
    #case (1)
      #local f_Lag = function(x,y,z) { 1 }
    #break
    #case (2)
      #switch (l)
        #case (0)
          #local n_Rad = n_Rad*2;
          #local f_Lag = function(x,y,z) { 1-rho(x,y,z) }
        #break
        #case (1)
          #local f_Lag = function(x,y,z) { 1 }
        #break
      #end
    #break
    #case (3)
      #switch (l)
        #case (0)
          #local f_Lag = function(x,y,z) { 3-6*rho(x,y,z)+2*pow(rho(x,y,z),2) }
        #break
        #case (1)
          #local n_Rad = n_Rad*2;
          #local f_Lag = function(x,y,z) { 2-rho(x,y,z) }
        #break
        #case (2)
          #local f_Lag = function(x,y,z) { 1 }
        #break
      #end
    #break
    #case (4)
      #switch (l)
        #case (0)
          #local n_Rad = n_Rad*4/3;
          #local f_Lag = function(x,y,z) { 3-9*rho(x,y,z)+6*pow(rho(x,y,z),2)-pow(rho(x,y,z),3) }
        #break
        #case (1)
          #local n_Rad = n_Rad*2;
          #local f_Lag = function(x,y,z) { 5-5*rho(x,y,z)+pow(rho(x,y,z),2) }
        #break
        #case (2)
          #local n_Rad = n_Rad*2;
          #local f_Lag = function(x,y,z) { 3-rho(x,y,z) }
        #break
        #case (3)
          #local f_Lag = function(x,y,z) { 1 }
        #break
      #end
    #break
  #end
  #declare f_Rad = function(x,y,z) { f_Lag(x,y,z)*exp(-rho(x,y,z)) }

  #declare Norm = n_Ang*n_Rad;

  #declare f_Orb = function(x,y,z) { f_Ang(x,y,z)*f_Rad(x,y,z) }
  #declare f_Den = function(x,y,z) { pow(f_Orb(x,y,z),2) }

  #declare Orb_pos = function(x,y,z) { value/Norm - f_Orb(x,y,z) }
  #declare Orb_neg = function(x,y,z) { value/Norm + f_Orb(x,y,z) }
  #declare Orb = function(x,y,z) { value/Norm - abs(f_Orb(x,y,z)) }
  #declare Den = function(x,y,z) { value/(Norm*Norm) - f_Den(x,y,z) }

#end