Quantum::ClebschGordan - Calculate/list Clebsch-Gordan Coefficients

    [davidrw@devbox davidrw]$ cg-j1j2 1 1/2

    use Quantum::ClebschGordan;

    print Quantum::ClebschGordan->new( j1=>1, j2=>'1/2', m=>'3/2', m1=>'1', m2=>'1/2', j=>'3/2' )->coeff;

    my $foo = Quantum::ClebschGordan->new();
    ...
    printf "%6s %6s %6s %6s %6s %6s %6s, %s\n", Quantum::ClebschGordan->state_names, 'N', 'c';
    printf "%6s %6s %6s %6s %6s %6s %6s, %s\n", $_->state_nums, $_->coeff, $_->coeff_value
      for Quantum::ClebschGordan->new( j1 => 1, j2 => '1/2' )->explode;

    perldoc Quantum::ClebschGordan

package Quantum::ClebschGordan;

use strict;
use warnings;

use base qw(Class::Accessor);

use Number::Fraction;
use Memoize;  # for factorial()
use Carp;

our $VERSION = '0.01';

our @STATE_VARS = qw/ j1 j2 m m1 m2 j /;
__PACKAGE__->mk_accessors( @STATE_VARS, '__coeff' );

sub new {
  my $self = shift;
  my $class = ref($self) || $self;
  my $p = { @_ };
  my $obj = bless {}, $class;
  foreach my $var ( @STATE_VARS ){
    $obj->set($var, $p->{$var} ) if exists $p->{$var};
  }
  return $obj;
}

sub set {
  my ($self, $key) = splice(@_, 0, 2);
  my @args = @_;
  my $settingStateVar = grep($key eq $_, @STATE_VARS) ? 1 : 0;
  if( $settingStateVar && scalar(@args) && defined $args[0] && ref($args[0]) ne 'Number::Fraction' ){
    my $v = Number::Fraction->new($args[0]);
    $args[0] = $v if defined $v;
  } 
  my $prev = $self->get($key);
  $self->SUPER::set($key, @args);
  if( $settingStateVar ){
    if( ! $self->__check_state ){
      $self->SUPER::set( $key, $prev );  # restore just in case this error is being trapped.
      croak "setting '$key' to " . join(",",@args) . " causes invalid state: " . join(',',@STATE_VARS) . " = " . join(',',map { defined $_ ? $_ : '' } $self->get(@STATE_VARS));
    }
    $self->__set_coeff();
  }
}

sub __check_state {
  # return undef on error
  # return 1 if all complete & valid
  # return -1 if incomplete but the parts given are valid
  my $self = shift;
  my ($j1, $j2, $m, $m1, $m2, $j) = $self->get(qw/j1 j2 m m1 m2 j/);
  carp("j1 is required") && return unless $j1;
  carp("j1 '$j1' is invalid (must be multiple of 1/2 and > 0)") && return unless $j1 > 0 && ($j1*2) == int ($j1*2);
  return -1 unless defined $j2;
  carp("j2 must be <= j1") && return unless $j2 <= $j1;
  carp("j2 '$j2' is invalid (must be multiple of 1/2 and > 0)") && return unless $j2 > 0 && $j2*2 == int $j2*2;
  return -1 unless defined $m;
  carp("m '$m' is not in range") && return unless -($j1+$j2) <= $m && $m <= ($j1+$j2);
  carp("m '$m' isn't valid in range") && return unless grep {$m == $_} map {-($j1+$j2) + $_} 0 .. 2*($j1+$j2);
  return -1 unless defined $m1;
  carp("m1 '$m1' is not in range") && return unless -$j1 <= $m1 && $m1 <= $j1;
  carp("m1 '$m1' isn't valid in range") && return unless grep {$m1 == $_} map {-$j1 + $_} 0 .. 2*$j1;
  return -1 unless defined $m2;
  carp("m2 '$m2' is not in range") && return unless -$j2 <= $m2 && $m2 <= $j2;
  carp("m2 '$m2' isn't valid in range") && return unless grep {$m2 == $_} map {-$j2 + $_} 0 .. 2*$j2;
  carp("m1+m2 != m for m=$m, m1=$m1, m2=$m2") && return unless $m == $m1 + $m2;
  return -1 unless defined $j;
  carp("j '$j' is not in range") && return unless abs($m) <= $j && $j <= $j1+$j2;
  carp("j '$j' isn't valid in range") && return unless grep {$j == $_} map {abs($m) + $_} 0 .. ($j1+$j2)-abs($m);
  return 1;
}

sub __seq {
  my ($start, $end, $iter) = @_;
  $iter = 1 if ! defined $iter;
  croak "'$iter' is not > 0" unless $iter && $iter > 0;
  $iter *= -1 if $start > $end;
  my $i = $start;
  my @list;
  push(@list, $i), $i+=$iter while ( ($iter>0) ? ($i <= $end) : ($i>= $end) );
  return @list;
}

# j1, j2
#  m = (j1+j2 down to -j1-j2)  by 1
#	2j+1 possibilities
#   m = m1 + m2
#   m1 = j1 .. -j1  by 1
#   m2 = j2 .. -j2  by 1
#  j = (j1+j2 down to abs(j1-j2) ) by 1
#  (2j1+1)(2j2+1) number of (j,m) pairs
#    j = j1+j2 down to abs(m)

sub explode {
  my $self = shift;
  croak "bad values" && return unless $self->__check_state;;
  my $j1 = $self->j1;
  my @j2 = defined $self->j2 ? ( $self->j2 ) : __seq(0.5, $j1, 0.5);
  my @coeffObjs;
  foreach my $j2 ( @j2 ){
    my @m = defined $self->m ? ( $self->m ) : __seq( $j1+$j2, -($j1+$j2) );
    my @m1 = defined $self->m1 ? ( $self->m1 ) : __seq( $j1, -$j1 );
    my @m2 = defined $self->m2 ? ( $self->m2 ) : __seq( $j2, -$j2 );
    foreach my $m ( @m ){
      my @j = defined $self->j ? ( $self->j ) : __seq( $j1+$j2, abs($m) );
      foreach my $m1 ( @m1 ){
        foreach my $m2 ( @m2 ){
          next unless $m == $m1 + $m2;
	  foreach my $j ( @j ){
#warn "j1 => $j1, j2 => $j2, m => $m, m1 => $m1, m2 => $m2, j => $j";
	    my $x = Quantum::ClebschGordan->new( j1 => $j1, j2 => $j2, m => $m, m1 => $m1, m2 => $m2, j => $j );
	    push @coeffObjs, $x;
#	    printf "%6s %6s %6s %6s %6s %6s %6s\n", $x->get(@STATE_VARS), $x->coeff;
	  }
        }
      }
    }
  }
  return @coeffObjs;
}

sub state_names {
  my $self = shift;
  return @STATE_VARS;
}

sub state_nums {
  my $self = shift;
  return $self->get(@STATE_VARS);
}

sub state_values {
  my $self = shift;
  return map { ref($_) eq 'Number::Fraction' ? $_->to_num : $_ } $self->get(@STATE_VARS);
}

memoize('factorial');
sub factorial {
  my $n = shift;
  return undef unless defined $n && $n =~ /^\d+$/; # non-negative integers only
  return 1 if $n == 0;
  return $n*factorial($n-1);
}


# From http://string.howard.edu/~tristan/QM2/QM2WE.pdf
#   "Addition of Angular Momenta, Clebsch-Gordan Coefficients and the Wigner-Eckart Theorem"
#   Section: '1.3. Clebsch-Gordan Coefficients'
#   Equations:
#     1.18a)	c(j,m, j1,j2,m1,m2) = delta(m,m1+m2) * rho(j, j1,j2) * sigma * tau
#     1.18b)	delta(m,m1+m2) =  (m==m1+m2) ? 1 : 0
#     1.18c)	rho = sqrt( (j1+j2-j)! * (j+j1-j2)! * (j2+j-j1)! * (2*j+1) / (j1+j2+j+1)! )
#     1.18d)	sigma = sqrt( (j+m)!*(j-m)!*(j1+m1)!*(j1-m1)!*(j2+m2)!*(j2-m2)! )
#     1.18e)	tau = SUM( -1^r / ( (j1-m1-r)! (j2+m2-r)! (j-j2+m1+r)! (j-j1-m2+r)! (j1+j2-j-r)! r! ) )
# 	note: 0! = 1; (-n)! = Gamma(1-n) = infinity for n = 1,2,...

sub __set_coeff {
  my $self = shift;
  $self->__coeff(undef);
  return unless $self->__check_state == 1;
  my ( $j1, $j2, $m, $m1, $m2, $j ) = $self->get( qw/ j1 j2 m m1 m2 j / );
  my $rho_squared_num = factorial($j1+$j2-$j) * factorial($j+$j1-$j2) * factorial($j2+$j-$j1) * (2*$j+1);
  my $rho_squared_den = factorial($j1+$j2+$j+1);
  my $sigma_squared = factorial($j+$m)*factorial($j-$m)*factorial($j1+$m1)*factorial($j1-$m1)*factorial($j2+$m2)*factorial($j2-$m2);

  # We'll store everything as a fraction to avoid floating point rounding errors.
  my $tau = Number::Fraction->new(0);
  # The sum over 4 is infinite, but for r's above this all the terms will be
  # zero (see note for equation 1.18e).
  my $r_max = 2*(abs($j)+abs($j1)+abs($j2)+abs($m1)+abs($m2));
  foreach my $r ( 0 .. $r_max ){
    my $denom;
    my @factorials = ( $j1-$m1-$r, $j2+$m2-$r, $j-$j2+$m1+$r, $j-$j1-$m2+$r, $j1+$j2-$j-$r, $r );
    foreach my $n ( @factorials ){
      my $n_factorial = factorial($n);
      undef $denom, last unless $n_factorial; # (-n)! is 1/0i (i.e. undef), so skip this term
      $denom = 1 unless defined $denom; # set to 1 the first time through.
      $denom *= $n_factorial;
    }
    next unless $denom;
    $tau += Number::Fraction->new( (-1)**$r , $denom );
  }
  # We want:
  # c(j,m, j1,j2,m1,m2) = delta(m,m1+m2) * rho(j, j1,j2) * sigma * tau
  #	= rho(j, j1,j2) * sigma * tau      # we've already bailed if delta=0
  #	= sign * sqrt( rho_squared * sigma_squared * tau * tau )
  #		where sign = sign(tau) ; we know rho & sigma are > 0
  #	= sign * sqrt( tau * tau * rho_squared * sigma_squared )  # flip to keep Number::Fraction objects
  #	= sign * sqrt( tau * tau * sigma_squared * rho_squared_num / rho_squared_den )
  my $sq = ($tau < 0 ? -1 : 1) * $tau * $tau * $sigma_squared * $rho_squared_num / $rho_squared_den;

  # still have a Number::Fraction object because of operator overloads.
  $self->__coeff( $sq->to_string );  # notation form
  
  return $self->coeff;
}

sub coeff {
  my $self = shift;
  return $self->__coeff;
}

sub coeff_value {
  my $self = shift;
  return notation2real( $self->__coeff );
}

sub notation2real {
  my $n = shift;
  $n or return $n; # covers 0, undef, ''
  $n = Number::Fraction->new( $n ) || $n;
  $n = $n->to_num if ref($n) eq 'Number::Fraction';
  my $sign = $n < 0 ? -1 : 1;
  return $sign * sqrt( abs($n) );
}

sub real2notation {
  my $x = shift;
  $x or return $x; # covers 0, undef, ''
  $x = $x->to_num if ref($x) eq 'Number::Fraction';
  my $sign = $x < 0 ? -1 : 1;
  my $multiplier = 10*9*8*7*6*5*4*3*2;
  $x = sprintf "%d", sprintf("%.4f",$x*$x*$multiplier);
  return $sign * Number::Fraction->new( $x, 1 ) / $multiplier;
}

sub wigner3j {
  my $self = shift;
  my $c = $self->coeff;
  return unless $self->__check_state() == 1;
  return $c unless $c;  # 0 or undef;
  $c = Number::Fraction->new($c) || $c;
  # (j1j2m1m2|j1j2m) = (-1)^(-j1+j2-m) * sqrt(2j+1) * Wigner(j1 j2 j, m1 m2 -m)
  my $power = -$self->j1 + $self->j2 - $self->m;
  $power = $power->to_num if ref($power) eq 'Number::Fraction';
  $c /= (-1)**( $power );
  my $sign = $c < 0 ? -1 : 1;
  $c /= 2*$self->j + 1;
  return $c;
}

sub wigner3j_value {
  my $self = shift;
  return notation2real( $self->wigner3j );
}

1; # End of Quantum::ClebschGordan