#include "global.h" //To compile : g++ AD.cxx -o {Input Executable Name} -lX11 //#include "Coeff.h" using namespace std; double RACAH(double A[6]){ double Racah = 0.; double ii[6],ix[6],kk[8]; double nn[4], nx[12],ny[8]; double K1,K2,K3,K4,K5,K6,K7,K8 = 0.; double NZMIN, NZMAX = .0; for(int i = 0;i<6;i++){ ii[i] = 2.0*A[i]; } double IA = ii[0]; double IB = ii[1]; double IC = ii[2]; double ID = ii[3]; double IE = ii[4]; double IF = ii[5]; K1 = IA + IB -IE; K3 = IC + ID -IE; K5 = IA + IC -IF; K7 = IB + ID -IF; K2 = IE - abs(IA-IB); K4 = IE - abs(IC-ID); K6 = IF - abs(IA-IC); K8 = IF - abs(IB-ID); double karr[8] = {K1,K2,K3,K4,K5,K6,K7,K8}; double kmn=karr[0]; double kmx=karr[0]; for(int i=0;i<8;i++) { if(kmn>karr[i]) { kmn=karr[i]; } else if(kmx<karr[i]) { kmx = karr[i]; } } double K0 = kmn; if(K0 < 0) return 0; /* 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); printf("K8 = %lf\n",K8); printf("K0 = %lf\n",K0); */ double nkk, n2kk,nkk2 = 0.; kk[0] = K1; kk[1] = K2; kk[2] = K3; kk[3] = K4; kk[4] = K5; kk[5] = K6; kk[6] = K7; kk[7] = K8; for(int i = 0; i < 7; i = i + 2){ nkk = kk[i]; n2kk = kk[i]/2; nkk2 = n2kk*2.; if(nkk != nkk2){return 0;} } ix[0] = IA; ix[1] = ID; ix[2] = IB; ix[3] = IC; ix[4] = IF; ix[5] = IE; double Lmin = ix[0]; int Kmin = 1; for( int i = 1; i < 6; i++){ if(ix[i] - Lmin < 0){ Lmin=ix[i]; } if(ix[i] - Lmin == 0){ Kmin=i; } if(ix[i] - Lmin > 0){ break; } } double sgn1 = 1.; double KQ,IP,IQ,IT = 0.; if(Kmin == 1){ KQ = 7 - Kmin; IA = ix[(int)KQ-1]; ID = ix[Kmin+4-1]; IE = ix[Kmin-1]; IF = ix[(int)KQ-4-1]; if(((int)(IA+ID-IE-IF)/2 %2) != 0) {sgn1 = -1.;} }else if(Kmin == 2){ KQ = 7 - Kmin; IA = ix[(int)KQ-1]; ID = ix[Kmin+4-1]; IE = ix[Kmin-1]; IF = ix[(int)KQ-4-1]; if(((int)(IA+ID-IE-IF)/2 %2) != 0) {sgn1 = -1.;} }else if(Kmin == 3){ KQ = 9 - Kmin; IB = ix[(int)KQ-1]; IC = ix[Kmin+2-1]; IE = ix[Kmin-1]; IF = ix[(int)KQ -2-1]; if(((int)(IB+IC-IE-IF)/2)%2 !=0) {sgn1 = -1.;} }else if(Kmin == 4){ KQ = 9 - Kmin; IB = ix[(int)KQ-1]; IC = ix[Kmin+2-1]; IE = ix[Kmin-1]; IF = ix[(int)KQ -2-1]; if(((int)(IB+IC-IE-IF)/2)%2 !=0) {sgn1 = -1.;} }else if(Kmin == 5){ IE = ix[4]; IF = ix[5]; IB = ix[3]; IC = ix[2]; }else if(Kmin == 6){ } IP = IA-IB; IQ = IC-ID; if(abs(IP) >= abs(IQ)){ if(IP >= 0){ IQ = IC-ID; double sgn2 = 1.0; if(((int)(IB+ID-IF)/2 % 2) != 0) { sgn2 = -1.;} if( IE > 0 ) { printf("Here 177\n"); 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); printf("it = %lf\n", IT); printf("ip = %lf\n", IP); printf("iq = %lf\n", IQ); // printf("bi = %lf\n", BI); nn[0] = (IA+IB+IE)/2 +1; nn[1] = (IC+ID+IE)/2 +1; nn[2] = (IA+IC+IF)/2 +1; nn[3] = (IB+ID+IF)/2 +1; nx[0] = nn[0] -IE; nx[1] = nn[0] -IB; nx[2] = nn[0] -IA; nx[3] = nn[1] -IE; nx[4] = nn[1] -ID; nx[5] = nn[1] -IC; nx[6] = nn[2] -IF; nx[7] = nn[2] -IC; nx[8] = nn[2] -IA; nx[9] = nn[3] -IF; nx[10] = nn[3] -ID; nx[11] = nn[3] -IB; IP = (IA+IB+IC+ID+4)/2; IQ = (IE+IF-IA-ID)/2; IT = (IE+IF-IB-IC)/2; double a = 1; double b = -IQ+1; if(a > b) { NZMIN = a; }else{ NZMIN = b; } double kzar[3] = {nx[1],nx[4],nx[10]}; double kzmn=kzar[0]; double kzmx=kzar[0]; for(int i=0;i<3;i++) { if(kzmn>kzar[i]) { kzmn=kzar[i]; } else if(kzmx<kzar[i]) { kzmx = kzar[i]; } } NZMAX = kzmn; printf("Here 237\n"); 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); printf("IT = %lf\n", IT); printf("IP = %lf\n", IP); printf("IQ = %lf\n", IQ); // printf("BI = %lf\n", BI); if(NZMAX < NZMIN){return 0;} double sqlog = 0; double I1 = 0.; for(int i = 0; i <4; i++){ I1 = nn[i]+1; sqlog = sqlog -lgamma(I1); } for(int i = 0; i <12; i++){ I1 = nx[i]; sqlog = sqlog - lgamma(I1); } sqlog = 0.5*sqlog; double sss = 0.; double NZM1 = 0.; double sslog = 0.; for(int nz = NZMIN; nz <= NZMAX; nz++){ NZM1 = nz - 1; ny[0] = IP - NZM1; ny[1] = nx[0] -NZM1; ny[2] = nx[3] -NZM1; ny[3] = nx[6] -NZM1; ny[4] = nx[9]-NZM1; ny[5] = nz; ny[6] = IQ+nz; ny[7] = IT+nz; I1 = ny[0]; sslog = sqlog + lgamma(I1); for(int i = 1; i<8; i++){ I1 = ny[i]; }sslog = sslog - lgamma(I1); sss = sss + (pow(-1.,NZM1)*exp(sslog)); } Racah = sgn1*sss; }else{ double BI = IB; double DI = ID; Racah = sgn1*sgn2/sqrt((BI+1.)*(DI+1.)); } }else{ IT = IB; IB = IA; IA = IT; IT = ID; ID = IC; IC = IT; IP = -IP; } }else{ IT = IC; IC = IA; IA = IT; IT = ID; ID = IB; IB = IT; IP = IQ; if(IP >= 0){ IQ = IC - ID; }else{ IT = IB; IB = IA; IA = IT; IT = ID; ID = IC; IC = IT; IP = -IP; } } double sgn2 = 1.; if((int)(IB+ID-IF)/2 %2 != 0) sgn2 = -1.; if(IE >= 0){ double BI = IB; double DI = ID; Racah = sgn1 * sgn2 / sqrt((BI+1.)*(DI+1.)); }else{ printf("Here 340\n"); 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); printf("IT = %lf\n", IT); printf("IP = %lf\n", IP); printf("IQ = %lf\n", IQ); // printf("BI = %lf\n", BI); nn[0] = (IA+IB+IE)/2 +1; nn[1] = (IC+ID+IE)/2 +1; nn[2] = (IA+IC+IF)/2 +1; nn[3] = (IB+ID+IF)/2 +1; nx[0] = nn[0] -IE; nx[1] = nn[0] -IB; nx[2] = nn[0] -IA; nx[3] = nn[1] -IE; nx[4] = nn[1] -ID; nx[5] = nn[1] -IC; nx[6] = nn[2] -IF; nx[7] = nn[2] -IC; nx[8] = nn[2] -IA; nx[9] = nn[3] -IF; nx[10] = nn[3] -ID; nx[11] = nn[3] -IB; IP = (IA+IB+IC+ID+4)/2; IQ = (IE+IF-IA-ID)/2; IT = (IE+IF-IB-IC)/2; double a = 1; double b = -IQ+1; if(a > b) { NZMIN = a; }else{ NZMIN = b; } double kzar[3] = {nx[1],nx[4],nx[10]}; double kzmn=kzar[0]; double kzmx=kzar[0]; for(int i=0;i<3;i++) { if(kzmn>kzar[i]) { kzmn=kzar[i]; } else if(kzmx<kzar[i]) { kzmx = kzar[i]; } } NZMAX = kzmn; printf("Here 401\n"); 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); printf("IT = %lf\n", IT); printf("IP = %lf\n", IP); printf("IQ = %lf\n", IQ); // printf("BI = %lf\n", BI); if(NZMAX < NZMIN){return 0;} double sqlog = 0; double I1 = 0.; for(int i = 0; i <4; i++){ I1 = nn[i]+1; sqlog = sqlog -lgamma(I1); } for(int i = 0; i <12; i++){ I1 = nx[i]; sqlog = sqlog - lgamma(I1); } sqlog = 0.5*sqlog; double sss = 0.; double NZM1 = 0.; double sslog = 0.; for(int nz = NZMIN; nz <= NZMAX; nz++){ NZM1 = nz - 1; ny[0] = IP - NZM1; ny[1] = nx[0] -NZM1; ny[2] = nx[3] -NZM1; ny[3] = nx[6] -NZM1; ny[4] = nx[9]-NZM1; ny[5] = nz; ny[6] = IQ+nz; ny[7] = IT+nz; I1 = ny[0]; sslog = sqlog + lgamma(I1); for(int i = 1; i<8; i++){ I1 = ny[i]; }sslog = sslog - lgamma(I1); sss = sss + (pow(-1.,NZM1)*exp(sslog)); } Racah = sgn1*sss; } return Racah; } int main(int argc, char** argv){ double j1 = 1.; double j2 = 2.; double l0abs = abs(j1-j2); double al0; if(l0abs > 0) { al0 = l0abs;} if(l0abs <= 0){ al0 = 1.;} double al1 = al0 + 1; double A[6] = {j1,j1,al0,al0,2,j2}; double B[6] = {j1,j1,al0,al0,4,j2}; double RCHA = RACAH(A); double RCHB = RACAH(B); printf("Racah A = %lf\n",RCHA); printf("Racah B = %lf\n",RCHB); return 0; }