815 lines
17 KiB
C++
815 lines
17 KiB
C++
|
|
#include "global.h"
|
|
//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
|
|
|
|
#include "Coeff.h"
|
|
|
|
using namespace std;
|
|
|
|
|
|
|
|
double CG2(double A[6]){
|
|
|
|
double ii[6];
|
|
double im[6];
|
|
double ix[9];
|
|
double nx[5];
|
|
|
|
for(int i = 0;i<6;i++){
|
|
ii[i] = 2.0*A[i];
|
|
}
|
|
double Cleb = 0.;
|
|
double IA = ii[0];
|
|
double IB = ii[1];
|
|
double IC = ii[2];
|
|
double ID = ii[3];
|
|
double IE = ii[4];
|
|
double IF = ii[5];
|
|
|
|
if(ID + IE - IF != 0){return 0;}
|
|
double K0 = IA + IB + IC;
|
|
if((int)K0 % 2 !=0){return 0;}
|
|
|
|
double K1 = IA + IB - IC;
|
|
double K2 = IC + IA - IB;
|
|
double K3 = IB + IC - IA;
|
|
double K4 = IA - abs(IB-IC);
|
|
double K5 = IB - abs(IC-IA);
|
|
double K6 = IC - abs(IA-IB);
|
|
|
|
double karr[6] = {K1,K2,K3,K4,K5,K6};
|
|
double kmn=karr[0];
|
|
double kmx=karr[0];
|
|
for(int i=0;i<6;i++)
|
|
{
|
|
if(kmn>karr[i])
|
|
{
|
|
kmn=karr[i];
|
|
}
|
|
else if(kmx<karr[i])
|
|
{
|
|
kmx = karr[i];
|
|
}
|
|
}
|
|
double K7 = kmn;
|
|
/*
|
|
printf("K1 = %lf\n",K1);
|
|
printf("K2 = %lf\n",K2);
|
|
printf("K3 = %lf\n",K3);
|
|
printf("K4 = %lf\n",K4);
|
|
printf("K5 = %lf\n",K5);
|
|
printf("K6 = %lf\n",K6);
|
|
printf("K7 = %lf\n",K7);
|
|
*/
|
|
|
|
if(K7 < 0) return 0;
|
|
|
|
for(int i = 0;i<3;i++){
|
|
if((int)(ii[i]+ii[i+3])%2 !=0)return 0;
|
|
if(ii[i] < abs(ii[i+3]))return 0;
|
|
}
|
|
double sgnfac =1.0;
|
|
|
|
for(int i = 0; i<6;i++){
|
|
im[i] = ii[i];
|
|
}
|
|
double FC2, IT = 0.;
|
|
//-------------------------------------------------------------------------
|
|
if(IA >= IB){
|
|
if(IC >= IB){
|
|
if(IB != 0){
|
|
if(IE < 0){
|
|
//assuming the logic passes;
|
|
FC2 = IC +1;
|
|
|
|
|
|
IT = K0/2 + 1;
|
|
ix[0] = IT-IC;
|
|
ix[1] = IT-IB;
|
|
ix[2] = IT-IA;
|
|
ix[3] = (IA+ID)/2 +1;
|
|
ix[4] = ix[3] - ID;
|
|
ix[5] = (IB+IE)/2 + 1;
|
|
ix[6] = ix[5] - IE;
|
|
ix[7] = (IC+IF)/2 +1;
|
|
ix[8] = ix[7] - IF;
|
|
//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
|
|
double sqfclg = log(FC2) - lgamma(IT+1);
|
|
|
|
|
|
|
|
double IXI;
|
|
for(int i = 0; i<9;i++){
|
|
IXI = ix[i];
|
|
sqfclg = sqfclg + lgamma(IXI);
|
|
}
|
|
sqfclg = 0.5 * sqfclg;
|
|
printf("Here place 1\n");
|
|
printf("ix[0] = %lf\n",ix[0]);
|
|
printf("ix[1] = %lf\n",ix[1]);
|
|
printf("ix[2] = %lf\n",ix[2]);
|
|
printf("ix[3] = %lf\n",ix[3]);
|
|
printf("ix[4] = %lf\n",ix[4]);
|
|
printf("ix[5] = %lf\n",ix[5]);
|
|
printf("ix[6] = %lf\n",ix[6]);
|
|
printf("ix[7] = %lf\n",ix[7]);
|
|
printf("ix[8] = %lf\n",ix[8]);
|
|
|
|
printf("sqdclg = %lf\n",sqfclg);
|
|
|
|
double xarr[3] = {ix[0],ix[4],ix[5]};
|
|
|
|
double xmn=xarr[0];
|
|
double xmx=xarr[0];
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
if(xmn>xarr[i])
|
|
{
|
|
xmn=xarr[i];
|
|
}
|
|
else if(xmx<xarr[i])
|
|
{
|
|
xmx = xarr[i];
|
|
}
|
|
}
|
|
double NZMAX = xmn;
|
|
|
|
double NZ2 = (IB-IC-ID)/2;
|
|
double NZ3 = (IA-IC+IE)/2;
|
|
|
|
double carr[3] = {0, NZ2, NZ3};
|
|
|
|
double cmn=carr[0];
|
|
double cmx=carr[0];
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
if(cmn>carr[i])
|
|
{
|
|
cmn=carr[i];
|
|
}
|
|
else if(cmx<carr[i])
|
|
{
|
|
cmx = carr[i];
|
|
}
|
|
}
|
|
|
|
|
|
double NZMIN = cmx+1;
|
|
|
|
if(NZMAX > NZMIN){return 0;}
|
|
double SS=0.;
|
|
double S1= pow((-1.),(NZMIN-1));
|
|
for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
|
|
|
|
double NZM1 = NZ -1;
|
|
nx[0] = ix[0]-NZM1;
|
|
nx[1] = ix[4]-NZM1;
|
|
nx[2] = ix[5]-NZM1;
|
|
nx[3] = NZ-NZ2;
|
|
nx[4] = NZ-NZ3;
|
|
double termlg = sqfclg - lgamma(NZ);
|
|
for(int i = 0; i<5; i++){
|
|
IXI = nx[i];
|
|
termlg= termlg -lgamma(IXI);
|
|
}
|
|
SS = SS + S1*exp(termlg);
|
|
S1 =-S1;
|
|
}
|
|
Cleb=sgnfac*SS;
|
|
|
|
|
|
}else{
|
|
// if (IE >=0) _>>>
|
|
ID = -ID;
|
|
IE = -IE;
|
|
IF = -IF;
|
|
if((int)(IA+IB-IC)/2 %2 != 0) sgnfac = (-1)*sgnfac;
|
|
|
|
IT = K0/2 + 1;
|
|
ix[0] = IT-IC;
|
|
ix[1] = IT-IB;
|
|
ix[2] = IT-IA;
|
|
ix[3] = (IA+ID)/2 +1;
|
|
ix[4] = ix[3] - ID;
|
|
ix[5] = (IB+IE)/2 + 1;
|
|
ix[6] = ix[5] - IE;
|
|
ix[7] = (IC+IF)/2 +1;
|
|
ix[8] = ix[7] - IF;
|
|
//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
|
|
FC2 = IC + 1;
|
|
double sqfclg = log(FC2) - lgamma(IT+1);
|
|
/*
|
|
printf("IT = %lf\n",IT);
|
|
printf("sqfclg = %lf\n",sqfclg);
|
|
|
|
printf("FC2 = %lf\n",FC2);
|
|
printf("lgamma(IT+1) = %lf\n",lgamma(IT+1));
|
|
printf("log(fC2) = %lf\n",log(FC2));
|
|
|
|
printf("IA = %lf\n",IA);
|
|
printf("IB = %lf\n",IB);
|
|
printf("IC = %lf\n",IC);
|
|
printf("ID = %lf\n",ID);
|
|
printf("IE = %lf\n",IE);
|
|
printf("IF = %lf\n",IF);
|
|
*/
|
|
double IXI;
|
|
for(int i = 0; i<9;i++){
|
|
IXI = ix[i];
|
|
sqfclg = sqfclg + lgamma(IXI);
|
|
// printf("sqfclg = %lf\n",sqfclg);
|
|
}
|
|
sqfclg = 0.5 * sqfclg;
|
|
/*
|
|
|
|
printf("Here place 2\n");
|
|
printf("ix[0] = %lf\n",ix[0]);
|
|
printf("ix[1] = %lf\n",ix[1]);
|
|
printf("ix[2] = %lf\n",ix[2]);
|
|
printf("ix[3] = %lf\n",ix[3]);
|
|
printf("ix[4] = %lf\n",ix[4]);
|
|
printf("ix[5] = %lf\n",ix[5]);
|
|
printf("ix[6] = %lf\n",ix[6]);
|
|
printf("ix[7] = %lf\n",ix[7]);
|
|
printf("ix[8] = %lf\n",ix[8]);
|
|
*/
|
|
printf("sqdclg = %lf\n",sqfclg);
|
|
double xarr[3] = {ix[0],ix[4],ix[5]};
|
|
|
|
double xmn=xarr[0];
|
|
double xmx=xarr[0];
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
if(xmn>xarr[i])
|
|
{
|
|
xmn=xarr[i];
|
|
}
|
|
else if(xmx<xarr[i])
|
|
{
|
|
xmx = xarr[i];
|
|
}
|
|
}
|
|
double NZMAX = xmn;
|
|
|
|
double NZ2 = (IB-IC-ID)/2;
|
|
double NZ3 = (IA-IC+IE)/2;
|
|
|
|
double carr[3] = {0, NZ2, NZ3};
|
|
|
|
double cmn=carr[0];
|
|
double cmx=carr[0];
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
if(cmn>carr[i])
|
|
{
|
|
cmn=carr[i];
|
|
}
|
|
else if(cmx<carr[i])
|
|
{
|
|
cmx = carr[i];
|
|
}
|
|
}
|
|
|
|
|
|
double NZMIN = cmx+1;
|
|
|
|
if(NZMAX > NZMIN){return 0;}
|
|
double SS=0.;
|
|
double S1= pow((-1.),(NZMIN-1));
|
|
|
|
printf("NZMAX = %lf\n",NZMAX);
|
|
printf("NZMIN = %lf\n",NZMIN);
|
|
printf("NZ2 = %lf\n",NZ2);
|
|
printf("NZ3 = %lf\n",NZ3);
|
|
printf("S1 = %lf\n",S1);
|
|
|
|
for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
|
|
|
|
double NZM1 = NZ -1;
|
|
nx[0] = ix[0]-NZM1;
|
|
nx[1] = ix[4]-NZM1;
|
|
nx[2] = ix[5]-NZM1;
|
|
nx[3] = NZ-NZ2;
|
|
nx[4] = NZ-NZ3;
|
|
double termlg = sqfclg - lgamma(NZ);
|
|
|
|
printf("nx[0] = %lf\n",nx[0]);
|
|
printf("nx[1] = %lf\n",nx[1]);
|
|
printf("nx[2] = %lf\n",nx[2]);
|
|
printf("nx[3] = %lf\n",nx[3]);
|
|
printf("nx[4] = %lf\n",nx[4]);
|
|
|
|
printf("termlg = %lf\n",termlg);
|
|
|
|
for(int i = 0; i<5; i++){
|
|
IXI = nx[i];
|
|
termlg= termlg -lgamma(IXI);
|
|
}
|
|
SS = SS + S1*exp(termlg);
|
|
S1 =-S1;
|
|
|
|
}
|
|
Cleb=sgnfac*SS;
|
|
|
|
}
|
|
// if (IB != 0) _>>>>
|
|
}else{
|
|
Cleb = sgnfac;
|
|
return Cleb;
|
|
|
|
}
|
|
// if (IC < IB) _>>>
|
|
|
|
}else{
|
|
|
|
IT = IC;
|
|
IC = IB;
|
|
IB = IT;
|
|
IT = IF;
|
|
IF = -IE;
|
|
IE = -IT;
|
|
sgnfac = sqrt((2.*A[2]+1.0)/(2.*A[1]+1.0));
|
|
if((int)(im[0]-im[3])/2 %2 != 0){sgnfac = (-1.)*sgnfac;}
|
|
|
|
Cleb = sgnfac;
|
|
return Cleb;
|
|
}
|
|
|
|
// if(IA < IB) _>>
|
|
}else{
|
|
|
|
if(IA >= IC){
|
|
IT = IC;
|
|
IC = IB;
|
|
IB = IT;
|
|
IT = IF;
|
|
IF = -IE;
|
|
IE = -IT;
|
|
|
|
Cleb = sgnfac;
|
|
return Cleb;
|
|
// if (IA < IC) as well as (IA < IB)
|
|
}else{
|
|
|
|
IT =IA;
|
|
IA =IB;
|
|
IB =IT;
|
|
IT =ID;
|
|
ID =IE;
|
|
IE =IT;
|
|
|
|
if((int)(K1/2) %2 != 0) {sgnfac = -1.;}
|
|
//call to 135 again.
|
|
if(IB != 0){
|
|
if(IE < 0){
|
|
//assuming the logic passes;
|
|
FC2 = IC +1;
|
|
|
|
|
|
IT = K0/2 + 1;
|
|
ix[0] = IT-IC;
|
|
ix[1] = IT-IB;
|
|
ix[2] = IT-IA;
|
|
ix[3] = (IA+ID)/2 +1;
|
|
ix[4] = ix[3] - ID;
|
|
ix[5] = (IB+IE)/2 + 1;
|
|
ix[6] = ix[5] - IE;
|
|
ix[7] = (IC+IF)/2 +1;
|
|
ix[8] = ix[7] - IF;
|
|
//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
|
|
double sqfclg = log(FC2) - lgamma(IT+1);
|
|
double IXI;
|
|
for(int i = 0; i<9;i++){
|
|
IXI = ix[i];
|
|
sqfclg = sqfclg + lgamma(IXI);
|
|
}
|
|
sqfclg = 0.5 * sqfclg;
|
|
printf("Here place 3\n");
|
|
printf("ix[0] = %lf\n",ix[0]);
|
|
printf("ix[1] = %lf\n",ix[1]);
|
|
printf("ix[2] = %lf\n",ix[2]);
|
|
printf("ix[3] = %lf\n",ix[3]);
|
|
printf("ix[4] = %lf\n",ix[4]);
|
|
printf("ix[5] = %lf\n",ix[5]);
|
|
printf("ix[6] = %lf\n",ix[6]);
|
|
printf("ix[7] = %lf\n",ix[7]);
|
|
printf("ix[8] = %lf\n",ix[8]);
|
|
|
|
printf("sqdclg = %lf\n",sqfclg);
|
|
double xarr[3] = {ix[0],ix[4],ix[5]};
|
|
|
|
double xmn=xarr[0];
|
|
double xmx=xarr[0];
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
if(xmn>xarr[i])
|
|
{
|
|
xmn=xarr[i];
|
|
}
|
|
else if(xmx<xarr[i])
|
|
{
|
|
xmx = xarr[i];
|
|
}
|
|
}
|
|
double NZMAX = xmn;
|
|
|
|
double NZ2 = (IB-IC-ID)/2;
|
|
double NZ3 = (IA-IC+IE)/2;
|
|
|
|
double carr[3] = {0, NZ2, NZ3};
|
|
|
|
double cmn=carr[0];
|
|
double cmx=carr[0];
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
if(cmn>carr[i])
|
|
{
|
|
cmn=carr[i];
|
|
}
|
|
else if(cmx<carr[i])
|
|
{
|
|
cmx = carr[i];
|
|
}
|
|
}
|
|
|
|
|
|
double NZMIN = cmx+1;
|
|
|
|
if(NZMAX > NZMIN){return 0;}
|
|
double SS=0.;
|
|
double S1= pow((-1.),(NZMIN-1));
|
|
for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
|
|
|
|
double NZM1 = NZ -1;
|
|
nx[0] = ix[0]-NZM1;
|
|
nx[1] = ix[4]-NZM1;
|
|
nx[2] = ix[5]-NZM1;
|
|
nx[3] = NZ-NZ2;
|
|
nx[4] = NZ-NZ3;
|
|
double termlg = sqfclg - lgamma(NZ);
|
|
for(int i = 0; i<5; i++){
|
|
IXI = nx[i];
|
|
termlg= termlg -lgamma(IXI);
|
|
}
|
|
SS = SS + S1*exp(termlg);
|
|
S1 =-S1;
|
|
|
|
}
|
|
Cleb=sgnfac*SS;
|
|
|
|
|
|
}else{
|
|
// if (IE >=0) _>>>
|
|
|
|
ID = -ID;
|
|
IE = -IE;
|
|
IF = -IF;
|
|
if((int)(IA+IB-IC)/2 %2 != 0) sgnfac = (-1.)*sgnfac;
|
|
|
|
IT = K0/2 + 1;
|
|
FC2 = IT +1;
|
|
ix[0] = IT-IC;
|
|
ix[1] = IT-IB;
|
|
ix[2] = IT-IA;
|
|
ix[3] = (IA+ID)/2 +1;
|
|
ix[4] = ix[3] - ID;
|
|
ix[5] = (IB+IE)/2 + 1;
|
|
ix[6] = ix[5] - IE;
|
|
ix[7] = (IC+IF)/2 +1;
|
|
ix[8] = ix[7] - IF;
|
|
//lgamma(arg) returns the log of the factorial of arg if arg is a natural number.
|
|
double sqfclg = log(FC2) - lgamma(IT+1);
|
|
double IXI;
|
|
for(int i = 0; i<9;i++){
|
|
IXI = ix[i];
|
|
sqfclg = sqfclg + lgamma(IXI);
|
|
}
|
|
sqfclg = 0.5 * sqfclg;
|
|
printf("Here place 4\n");
|
|
printf("ix[0] = %lf\n",ix[0]);
|
|
printf("ix[1] = %lf\n",ix[1]);
|
|
printf("ix[2] = %lf\n",ix[2]);
|
|
printf("ix[3] = %lf\n",ix[3]);
|
|
printf("ix[4] = %lf\n",ix[4]);
|
|
printf("ix[5] = %lf\n",ix[5]);
|
|
printf("ix[6] = %lf\n",ix[6]);
|
|
printf("ix[7] = %lf\n",ix[7]);
|
|
printf("ix[8] = %lf\n",ix[8]);
|
|
|
|
printf("sqdclg = %lf\n",sqfclg);
|
|
double xarr[3] = {ix[0],ix[4],ix[5]};
|
|
|
|
double xmn=xarr[0];
|
|
double xmx=xarr[0];
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
if(xmn>xarr[i])
|
|
{
|
|
xmn=xarr[i];
|
|
}
|
|
else if(xmx<xarr[i])
|
|
{
|
|
xmx = xarr[i];
|
|
}
|
|
}
|
|
double NZMAX = xmn;
|
|
|
|
double NZ2 = (IB-IC-ID)/2;
|
|
double NZ3 = (IA-IC+IE)/2;
|
|
|
|
double carr[3] = {0, NZ2, NZ3};
|
|
|
|
double cmn=carr[0];
|
|
double cmx=carr[0];
|
|
for(int i=0;i<3;i++)
|
|
{
|
|
if(cmn>carr[i])
|
|
{
|
|
cmn=carr[i];
|
|
}
|
|
else if(cmx<carr[i])
|
|
{
|
|
cmx = carr[i];
|
|
}
|
|
}
|
|
|
|
|
|
double NZMIN = cmx+1;
|
|
|
|
if(NZMAX > NZMIN){return 0;}
|
|
double SS=0.;
|
|
double S1= pow((-1.),(NZMIN-1));
|
|
for(int NZ = NZMIN;NZ<=NZMAX;NZ++){
|
|
|
|
double NZM1 = NZ -1;
|
|
nx[0] = ix[0]-NZM1;
|
|
nx[1] = ix[4]-NZM1;
|
|
nx[2] = ix[5]-NZM1;
|
|
nx[3] = NZ-NZ2;
|
|
nx[4] = NZ-NZ3;
|
|
double termlg = sqfclg - lgamma(NZ);
|
|
for(int i = 0; i<5; i++){
|
|
IXI = nx[i];
|
|
termlg= termlg -lgamma(IXI);
|
|
}
|
|
SS = SS + S1*exp(termlg);
|
|
S1 =-S1;
|
|
}
|
|
Cleb=sgnfac*SS;
|
|
|
|
}
|
|
// if (IB != 0) _>>>>
|
|
}else{
|
|
Cleb = sgnfac;
|
|
return Cleb;
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
return Cleb;
|
|
}
|
|
|
|
|
|
//double CG(double J, double M, double j1, double m1, double j2, double m2){
|
|
double CG(double j1, double j2, double J, double m1, double m2, double M){
|
|
//recall that j1,m1 + j2,m2 = J,M
|
|
|
|
printf("-----------------\n");
|
|
if(M != m1 + m2) return 0;
|
|
|
|
double Jmin = abs(j1 - j2);
|
|
double Jmax = j1+j2;
|
|
|
|
printf(" Jmin = %lf\n", Jmin);
|
|
printf(" Jmax = %lf\n", Jmax);
|
|
|
|
if(J < Jmin || Jmax < J) return 0;
|
|
|
|
// double a0 = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0);
|
|
|
|
|
|
double a0 = (2*J+1.0)*tgamma(J+j1-j2+1) * tgamma(J-j1+j2+1) * tgamma(j1+j2-J+1)/tgamma(J+j1+j2+1.0 +1);
|
|
double A0 = sqrt(a0);
|
|
|
|
// double a = factorial(J+M) *factorial(J-M);
|
|
// double a1= factorial(j1+m1) *factorial(j1-m1);
|
|
// double a2= factorial(j2+m2) *factorial(j2-m2);
|
|
|
|
double a = tgamma(J+M+1) *tgamma(J-M+1);
|
|
double a1= tgamma(j1+m1+1) *tgamma(j1-m1+1);
|
|
double a2= tgamma(j2+m2+1) *tgamma(j2-m2+1);
|
|
|
|
double A = sqrt( a * a1 * a2);
|
|
|
|
printf(" a0 = %lf\n", a0);
|
|
printf(" A0 = %lf\n", A0);
|
|
printf(" a = %lf\n", a);
|
|
printf(" a1 = %lf\n", a1);
|
|
printf(" a2 = %lf\n", a2);
|
|
printf(" A = %lf\n", A);
|
|
|
|
int pmax = min( min(j1+j2-J,j1-m1),j2 + m2);
|
|
|
|
printf("pmax = %d\n",pmax);
|
|
|
|
double cg = 0.;
|
|
|
|
for( int p =0; p<=pmax;p++){
|
|
|
|
// double p1 = factorial(j1+j2-J-p);
|
|
// double p2 = factorial(j1-m1-p);
|
|
// double p3 = factorial(j2+m2-p);
|
|
// double p4 = factorial(J -j2 +m1 +p);
|
|
// double p5 = factorial(J -j1 -m2 +p);
|
|
// double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5);
|
|
|
|
|
|
double p1 = tgamma(j1+j2-J-p+1);
|
|
double p2 = tgamma(j1-m1-p+1);
|
|
double p3 = tgamma(j2+m2-p+1);
|
|
double p4 = tgamma(J -j2 +m1 +p+1);
|
|
double p5 = tgamma(J -j1 -m2 +p+1);
|
|
double t = pow(-1,p)/(tgamma(p+1) * p1 * p2 * p3 * p4 * p5);
|
|
|
|
printf("t = %lf\n",t);
|
|
cg += t;
|
|
printf("cg = %lf\n",cg);
|
|
}
|
|
printf("-----------------\n");
|
|
return A0*A*cg;
|
|
}
|
|
|
|
double CG_a(double A[6]){
|
|
//recall that j1,m1 + j2,m2 = J,M
|
|
//testing assingments until the output is the same.
|
|
double j1 = A[0];
|
|
double j2 = A[1];
|
|
double J = A[2];
|
|
double m1 = A[3];
|
|
double m2 = A[4];
|
|
double M = A[5];
|
|
|
|
|
|
|
|
if(M != m1 + m2) return 0;
|
|
|
|
double Jmin = abs(j1 - j2);
|
|
double Jmax = j1+j2;
|
|
|
|
if(J < Jmin || Jmax < J) return 0;
|
|
|
|
double a0 = (2*J+1.0)*factorial(J+j1-j2) * factorial(J-j1+j2) * factorial(j1+j2-J)/factorial(J+j1+j2+1.0);
|
|
double A0 = sqrt(a0);
|
|
|
|
double a = factorial(J+M) *factorial(J-M);
|
|
double a1= factorial(j1+m1) *factorial(j1-m1);
|
|
double a2= factorial(j2+m2) *factorial(j2-m2);
|
|
double A1 = sqrt( a * a1 * a2);
|
|
|
|
int pmax = min( min(j1+j2-J,j1-m1),j2 + m2);
|
|
|
|
double cg = 0.;
|
|
|
|
for( int p =0; p<=pmax;p++){
|
|
|
|
double p1 = factorial(j1+j2-J-p);
|
|
double p2 = factorial(j1-m1-p);
|
|
double p3 = factorial(j2+m2-p);
|
|
double p4 = factorial(J -j2 +m1 +p);
|
|
double p5 = factorial(J -j1 -m2 +p);
|
|
double t = pow(-1,p)/(factorial(p) * p1 * p2 * p3 * p4 * p5);
|
|
|
|
cg += t;
|
|
}
|
|
return A0*A1*cg;
|
|
}
|
|
double ThreeJsym(double j1, double m1, double j2, double m2, double j3, double m3){
|
|
//[j1 j2 j3] = (-1)^(j1-j2-m3)/ sqrt(2*j3+1) * CG(j3, -m3, j1, m1, j2, m2)
|
|
//[m1 m2 m3]
|
|
|
|
|
|
|
|
return pow(-1,j1 -j2 -m3)/sqrt(2*j3+1)*CG(j3,-m3,j1,m1,j2,m2);
|
|
}
|
|
|
|
double SixJsym(double j1, double j2, double j3, double j4, double j5, double j6){
|
|
|
|
//--------------------------------------------------------------------------------//
|
|
// The six j symbol describes the coupling between j1 j2 and j3 to J - j1.
|
|
// essentially a triangle of angular momentum rules between these.
|
|
// j1 = j1
|
|
// j2 = j2
|
|
// j3 = j1 + j2
|
|
// j4 = j3
|
|
// j5 = J = j1 + j2 + j3
|
|
// j6 = j2 + j3
|
|
// ----------------------------------------------------------------------------- //
|
|
// the following conditions check the triangle selection rules
|
|
|
|
if( j3 < abs(j1 - j2) || j1 + j2 < j3) return 0;
|
|
if( j6 < abs(j2 - j4) || j2 + j4 < j6) return 0;
|
|
if( j5 < abs(j1 - j6) || j1 + j6 < j5) return 0;
|
|
if( j5 < abs(j3 - j4) || j3 + j4 < j5) return 0;
|
|
// now that they have been checked, we can go ahead and calculate sixJ.
|
|
|
|
double sixj = 0.0;
|
|
float m1 = -j1;
|
|
float m2 = -j2;
|
|
float m3 = -j3;
|
|
float m4 = -j4;
|
|
float m5 = -j5;
|
|
float m6 = -j6;
|
|
|
|
for(; m1 <= j1; m1 = m1 +1){
|
|
for(; m2 <= j2; m2 = m2 +1){
|
|
for(; m3 <= j3; m3 = m3 + 1){
|
|
for(; m4 <= j4; m4 = m4 +1){
|
|
for(; m5 <= j5; m5 = m5 + 1){
|
|
for(; m6 <= j6; m6 = m6 +1){
|
|
|
|
double h = (j1 - m1) + (j2 - m2) + (j3 -m3) + (j4 - m4) + (j5 - m5) + (j6 - m6);
|
|
double b1 = ThreeJsym(j1, -m1, j2, -m2, j3, -m3);
|
|
double b2 = ThreeJsym(j1, m1, j5, -m5, j6, m6);
|
|
double b3 = ThreeJsym(j4, m4, j2, m2, j6, -m6);
|
|
double b4 = ThreeJsym(j4, -m4, j5, m5, j3, m3);
|
|
double b = b1 * b2 * b3 * b4;
|
|
|
|
sixj += pow(-1,h)*b;
|
|
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return sixj;
|
|
}
|
|
|
|
|
|
int main(int argc, char ** argv){
|
|
|
|
//if mod (2*J1,2) = 1 do this
|
|
// double A0[6] = {1,1,2,0,0,0};
|
|
// double A1[6] = {1,1,4,0,0,0};
|
|
|
|
//if mod(2*j1,20 = 0 do this
|
|
|
|
double j1 = 1.;
|
|
double j2 = 2.;
|
|
|
|
|
|
double A0[6] = {j1,j1,2,.5,-.5,0};
|
|
double A1[6] = {j1,j1,2,-.5,.5,0};
|
|
double A2[6] = {j1,j1,4,-.5,.5,0};
|
|
double A3[6] = {j1,j1,4,.5,-.5,0};
|
|
|
|
/*
|
|
double cg0 = CG2(A0);
|
|
double cg1 = CG2(A1);
|
|
double cg2 = CG2(A2);
|
|
double cg3 = CG2(A3);
|
|
*/
|
|
double cgr0= CG(j1,j1,2,.5,-.5,0);
|
|
double cgr1= CG(j1,j1,2,-.5,.5,0);
|
|
double cgr2= CG(j1,j1,4,.5,-.5,0);
|
|
double cgr3= CG(j1,j1,4,-.5,.5,0);
|
|
/*
|
|
printf("------\n");
|
|
printf("CG0 = %lf\n",cg0);
|
|
printf("CG1 = %lf\n",cg1);
|
|
printf("CG2 = %lf\n",cg2);
|
|
printf("CG3 = %lf\n",cg3);
|
|
*/
|
|
printf("------\n");
|
|
printf("CGR0 = %lf\n",cgr0);
|
|
printf("CGR1 = %lf\n",cgr1);
|
|
printf("CGR2 = %lf\n",cgr2);
|
|
printf("CGR3 = %lf\n",cgr3);
|
|
|
|
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
|