diff --git a/1062.csv b/1062.csv
new file mode 100644
index 0000000..98dd2da
--- /dev/null
+++ b/1062.csv
@@ -0,0 +1,4 @@
+angle,Y,Yerr
+45,129,10
+90,110,10
+135,129,10
diff --git a/AD.cxx b/AD.cxx
new file mode 100644
index 0000000..368144d
--- /dev/null
+++ b/AD.cxx
@@ -0,0 +1,750 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+
+#include "GUI.h"
+#include "Coeff.h"
+#include "QDK.h"
+#include "Racah.h"
+#include "Cleb.h"
+
+
+//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
+//For more information about angular distributions, read Rose and Brinks, and Frank Moore (1988). 
+
+using namespace std;
+void menu(){
+
+	std::cout <<"==========================================\n";
+	std::cout<< "\t\t Menu Options \t \n";
+	std::cout<< "0 - Reading and Instructions \n"; 
+	std::cout<< "1 - Plot Chi Sqr \n";
+	std::cout<< "2 - Plot Ang.Dis Fit \n"; 
+	std::cout<< "3 - Clear J1 J2 Memory \n";
+	std::cout<< "4 - Exit \n";
+	std::cout<< "5 - Legfit to Data \n"; 
+	std::cout<< "6 - Generate Data \n"; 
+	std::cout<<"=========================================== \n";} 
+void Readme(){
+	std::cout<< " The program calculates Chi-Squared values \n";
+	std::cout<< " from experimental angular distributions as \n";
+	std::cout<< " a function of multipole ratios using the \n"; 
+	std::cout<< " theoretical angular distribution formulae \n"; 
+	std::cout<< " in rose and brink \n";
+	std::cout<< " \n"; 
+	std::cout<< " Follow the prompt in order to correctly display \n";
+	std::cout<< " \n";
+	std::cout<< " To close the gui, press any button. \n";
+	std::cout<< " To zoom in, left click then drag and let go. To unzoom\n";
+	std::cout<< " press the space bar. To draw, right click\n";
+	std::cout<< "                                          \n";
+	std::cout<< " ad.txt is generated with geometric stats.\n";
+	std::cout<< " Then at the end it generates CG, W3J & W6J\n";
+	
+}
+#ifndef __CINT__ 
+
+int main(int argc,char **argv){
+	HistoGUI gui;
+	//	menu();
+	double j1, j2;
+
+
+	//detector data input
+	double detradius, targetdistance, detthickness = -1;
+	double Energy;
+	//input .csv file
+	double Sigma; 	
+	double Feeding; 
+
+//--------------------------------
+
+	detradius = 3.;
+	targetdistance = 4.;
+	detthickness = 5.;
+	
+	Energy = 1062.;
+	Sigma = .1; 
+	Feeding = .1;
+
+	j1 = 2.; 
+	j2 = 1.;
+//--------------------------------
+	
+//TEST MODE? y : 1 || n = 0
+
+	int test = 1; 
+
+	int det_param_token = -1;
+	int gamma_energy_token = -1;
+	int ang_file_token = -1; 
+	int sigma_token = -1; 
+	int feeding_token = -1;
+	int j1j2token = -1; 
+ if(test == 1){
+
+	 det_param_token = 1;
+	 gamma_energy_token = 1;
+	 ang_file_token = -1; 
+	 sigma_token = 1; 
+	 feeding_token = 1;
+	 j1j2token = 1; 
+
+}
+/*
+	printf("Input detector radius, target distance, and detector thickness : ");
+	scanf("%lf,%lf,%lf",&detradius,&targetdistance,&detthickness);
+
+
+	if(detradius > 0 && targetdistance > 0 && detthickness > 0){
+		det_param_token = 1;
+	}else{
+
+		do{
+			if(detradius< 0|| targetdistance < 0|| detthickness < 0){
+				printf("Negative numbers are not allowed!\nRe-enter radius, distance, and thickness : ");
+				scanf("%lf,%lf,%lf", &detradius,&targetdistance,&detthickness);
+			}
+		}while(detradius< 0|| targetdistance < 0|| detthickness < 0);}
+
+	if(detradius > 0 && targetdistance > 0 && detthickness > 0){
+		det_param_token = 1; printf(" --- Detector characteristics loaded --- \n");
+	}
+
+	printf("Enter the gamma-ray energy : ");
+	scanf("%lf", &Energy);
+	if(Energy > 0){
+		gamma_energy_token = 1;
+	}else{
+		do{
+			if(Energy <= 0){		
+				printf("Negative numbers are not allowed!\nRe-enter gamma-ray energy : ");
+			scanf("%lf", &Energy);}
+
+		}while(Energy <= 0);}
+
+	if(Energy > 0){
+		gamma_energy_token = 1; printf(" --- Gamma Energy loaded --- \n");
+	}*/
+
+	QKN(Energy, detradius, targetdistance, detthickness);
+	//then calucate QD2 and QD4, and replace the 0 with them. 
+	double QD2 = .7;
+	double QD4 = .3;
+
+	//input data file of Angular Data (theta, Yexp, Yerr)
+	string fname;
+	cout<<"Enter the file name : ";
+	cin>>fname;
+
+	vector<vector<string> > content;
+	vector<string> row;
+	string line, word; 
+	
+	vector<string> adata;
+	vector<string> ydata;
+	vector<string> eydata; 
+	int num = 0;
+	fstream file(fname.c_str());
+	if(file.is_open()){
+		while(getline(file,line)){
+			row.clear();
+	
+			stringstream str(line);
+
+			while(getline(str, word, ','))
+			
+				row.push_back(word);
+				content.push_back(row);
+				adata.push_back(row[0]);
+				ydata.push_back(row[1]);
+				eydata.push_back(row[2]);			
+				 
+		}
+		cout<< " --- Angular Data File Loaded --- \n";
+	}else cout<< "Could not open the file\n";
+	
+	for(int i = 0; i< adata.size();i++){
+//		cout << adata[i] << "\n";
+//		cout << ydata[i] << "\n";
+//		cout << eydata[i] << "\n";
+	}
+	
+	// Legendre Polinomial fit using adata, ydata, eydata;
+
+	vector<int> int_adata(adata.size());
+std::transform(adata.begin()+1, adata.end(), std::back_inserter(int_adata), StrToInt);
+
+	vector<double> double_angle(int_adata.size() -4);
+	for(int i = 4; i<  int_adata.size(); i++){
+	//	printf("angle = %d\n",int_adata[i]);
+		double_angle.push_back((double)int_adata[i]*3.14159/180);	
+	}
+
+	vector<int> int_ydata(ydata.size());
+std::transform(ydata.begin()+1, ydata.end(), std::back_inserter(int_ydata), StrToInt);
+
+	for(int i = 4; i<  int_ydata.size(); i++){
+	//	printf("Ydata = %d\n",int_ydata[i]);
+	}
+
+	vector<double> double_ydata(int_ydata.size() -4);
+	for(int i = 4; i<  int_ydata.size(); i++){
+	//	printf("angle = %d\n",int_adata[i]);
+		double_ydata.push_back((double)int_ydata[i]);	
+
+	}
+	
+	vector<int> int_eydata(eydata.size());
+std::transform(eydata.begin()+1, eydata.end(), std::back_inserter(int_eydata), StrToInt);
+
+	for(int i = 4; i<  int_eydata.size(); i++){
+	//	printf("Error data = %d\n",int_eydata[i]);
+	}
+
+	
+	for(int i = 3; i< double_angle.size();i++){
+	//	printf("rad angle = %lf\n",double_angle[i]);
+	
+		//displays the angles inputed as doubles in radians. 
+	}
+	
+	for(int i = 3; i< double_ydata.size();i++){
+	//	printf("y_data = %lf\n",double_ydata[i]);
+	
+		//displays the angles inputed as doubles in radians. 
+	}
+	/*
+	for(int i=0; i<content.size(); i++){
+		for(int j=0; j<content[i].size(); j++){
+			cout<< content[i][j]<<" ";
+		}
+		cout<<"\n";
+	}
+	*/
+	//legendre fitting 
+	// F[xi] = A0 + A2(P2(xi)) + A4(P4(xi))
+	//xi data is angle data in radians. 
+	//yi data is y-intensity data. 
+	
+	//to change the number of angles, just mess with the indexing of the read in and arrays. 
+
+
+	//12 constants. this will build our matrix for gaussian elinmination. 
+	//[a1 a2   a3] [A0] = [a4]
+	//[a5 a6   a7] [A2] = [a8]
+	//[a9 a10 a11] [A4] = [a12]
+	
+	double a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12 = 0.;
+	
+	for(int i = 3; i< double_angle.size(); i++){
+		
+//eq1
+		a1 += 1; //how many data points, not intensional but convientient.  
+		a2 += (1.5 * pow(double_angle[i],2) - 1);
+		a3 += (35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a4 += double_ydata[i]; 
+//eq2 
+		a5 += (1.5 * pow(double_angle[i],2) - 1);
+		a6 += (1.5 * pow(double_angle[i],2) - 1)*(1.5 * pow(double_angle[i],2) - 1);
+		a7 += (1.5 * pow(double_angle[i],2) - 1)*(35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a8 += (1.5 * pow(double_angle[i],2) - 1)*double_ydata[i];
+//eq4
+		a9 += (35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a10 += (1.5 * pow(double_angle[i],2) - 1)*(35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a11 +=(35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. )* (35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. );
+		a12 += (35./8. * pow(double_angle[i],4) - 30./8. * pow(double_angle[i],2)  +  3./8. )*double_ydata[i];
+
+		 
+	}
+/*
+	printf("a1 = %lf\n",a1);	
+	printf("a2 = %lf\n",a2);	
+	printf("a3 = %lf\n",a3);	
+	printf("a4 = %lf\n",a4);	
+	printf("a5 = %lf\n",a5);	
+	printf("a6 = %lf\n",a6);	
+	printf("a7 = %lf\n",a7);	
+	printf("a8 = %lf\n",a8);	
+	printf("a9 = %lf\n",a9);	
+	printf("a10 = %lf\n",a10);	
+	printf("a11 = %lf\n",a11);	
+	printf("a12 = %lf\n",a12);	
+*/
+//here number of equations are 3 because its a 3 coeff. fit. 
+	int n = 3; 
+	double mat[3][4] = {{a1,a2,a3,a4},{a5,a6,a7,a8},{a9,a10,a11,a12}};
+//where A0,A2,A4 are stored
+	double residual[n];
+
+	for(int i = 0; i<n;i++){
+		for( int j = i+1; j<n; j++){
+			if(abs(mat[i][i]) < abs(mat[i][j])){
+				for(int k = 0; k<n+1;k++){
+					//swapping matrix
+		//			mat[i][k] = mat[i][k] + mat[j][k];
+		//			mat[j][k] = mat[i][k] - mat[j][k];
+		//			mat[i][k] = mat[i][k] - mat[j][k];
+					swap(mat[i][k],mat[j][k]);
+				}
+
+			}
+		}
+	}
+//performing gaussian elimination.
+
+	for(int i = 0; i< n-1; i++) {
+
+		for(int j = i +1; j<n;j++){
+			
+			float f = mat[j][i]/mat[i][i];
+
+			for(int k = 0; k<n+1; k++){
+				mat[j][k] = mat[j][k] -f*mat[i][k];
+
+			}
+		}
+	}
+	
+//backwards substitution for unknowns.
+
+	for( int i = n-1; i >= 0; i--){
+
+		residual[i] = mat[i][n];
+
+		for( int j = i + 1; j<n;j++){
+
+			if(i != j){
+				residual[i] = residual[i] - mat[i][j] * residual[j];
+			}
+		}
+		residual[i] = residual[i]/mat[i][i]; 
+	}
+	
+
+	cout<<"A0 = "<<residual[0]<<"\n";
+	cout<<"A2 = "<<residual[1]<<"\n";
+	cout<<"A4 = "<<residual[2]<<"\n";
+
+	
+
+
+	/*
+
+	printf("Please enter sigma, sigma = 0 for perfect allignment : ");
+	scanf("%lf", &Sigma);
+	if(Sigma >= 0 ){
+		sigma_token = 1; 
+	}else{
+		do{
+
+			printf("Negative numbers are not allowed!\nRe-enter Sigma : ");
+			scanf("%lf", &Sigma);
+		}while(Sigma <0);}
+
+	if(Sigma >= 0 ){
+		sigma_token = 1; printf(" --- Sigma Loaded --- \n"); 
+	}
+
+	printf("Please enter feeding parameter : ");
+	scanf("%lf", &Feeding);	
+	if(Feeding >= 0 and Feeding < 1){
+		feeding_token = 1;
+	}else{
+		do{
+			printf("Negative numbers are not allowed!\nRe-enter Feeding : ");
+			scanf("%lf", &Feeding);
+		}while(Feeding < 0);}
+
+	if(Feeding >= 0 ){
+		feeding_token = 1; printf(" --- Feeding Loaded --- \n");
+	}
+
+	printf("Please enter J1,J2 : ");
+	scanf("%lf,%lf",&j1,&j2);	
+	if(j1 >= 0 && j2  >= 0 ){
+		j1j2token = 1;
+	}else{
+		do{
+			printf("Negative numbers are not allowed!\nRe-enter J1,J2 : ");
+			scanf("%lf,%lf", &j1,&j2);	
+		}while(j1 < 0 || j2 < 0);}
+
+	if(j1 >= 0 && j2  >= 0 ){
+		j1j2token = 1; printf(" --- J1 & J2 Loaded --- \n");
+	}
+*/
+	
+
+	/*
+	   printf("detector radius = %lf\n",detradius);
+	   printf("detector distance = %lf\n",targetdistance);
+	   printf("detector width = %lf\n",detthickness);
+
+	   printf("Gamma-Ray energy = %lf\n", Energy);
+	   printf("Sigma Distribution = %lf\n", Sigma); 
+	   printf("Feeding param = %lf\n", Feeding); 
+	   printf("J1,J2 = %lf,%lf\n",j1,j2);
+	   */
+/*	
+	double J  = j1 + j2;
+	double m1 = j1;
+	double m2 = j2;
+	double M  = m1 + m2; 
+	double j3 = J; 
+	double m3 = j3;	
+	double j4 = j3;
+	double j5 = j1 + j2 + j3;
+	double j6 = j2 + j3; 
+*/
+// Calculate BK(j1) for perfect allignment
+//
+	double A[6];
+	double js = sqrt(2*j1 +1); 
+	double j12 = 2*j1;
+	int IS = (int)j12 % 2; 
+	
+	double cg12,cg22,cg24,cg14 = 0.;
+	double I1, I2 = 0.;
+	double Bk11, Bk12 = 0.;
+	
+	double cg1, cg2 = 0.;	
+	
+	if(IS == 1){
+		A[0] = j1;
+		A[1] = j1;
+		A[2] = 2.;
+		A[3] = .5;
+		A[4] = -.5;
+		A[5] = 0.;
+	
+		cg12 = CG2(A);
+		
+		A[3] = -.5;
+		A[4] = .5;
+		
+		cg22 = CG2(A);
+
+		A[2] = 4.;
+		
+		cg24 = CG2(A);
+
+		A[3] = .5;
+		A[4] = -.5;
+
+		cg14 = CG2(A);
+
+		Bk11 = (0.5)*js*(pow(-1,I1) * cg12 + pow(-1,I2) * cg22);
+		Bk12 = (0.5)*js*(pow(-1,I1) * cg14 + pow(-1,I2) * cg24);
+	}else if(IS == 0){
+		A[0] = j1; 
+		A[1] = j1; 
+		A[2] = 2.;
+		A[3] = 0.;
+		A[4] = 0.;
+		A[5] = 0.; 
+		
+		cg1 = CG2(A);
+	
+		Bk11 = pow(-1,j1) * js * cg1;
+
+		A[2] = 4.;
+	
+		cg2 = CG2(A);
+		Bk12 = pow(-1,j1) * js * cg2; 	
+		
+	}
+
+//normalize W(m1) 
+//
+
+	       j12 = 2*j1;
+	double j14 = 4*j1;
+	
+	double sigsq = pow(Sigma,2);
+	double sum1 = 0.; 
+
+	double am1,amsq,x,ex = 0.;
+
+	for(int i = 0; i <= j14; i = i + 2){
+		am1 = 0.5*(i - j12);
+		amsq = pow(am1,2);
+		x = - (amsq/(2*sigsq));
+		ex = exp(x);
+
+		sum1 = sum1 + ex;
+	}
+	double cn1 = 1./sum1;	
+
+	double AL0 = 0.;
+	double A0 = j1 - j2; 
+	double A0b = abs(A0);
+	if(A0b > 0){ AL0 = A0b;}
+	
+	if(A0b <= 0){AL0 = 1.;}
+
+	double AL1 = AL0 +1; 
+
+	double am11,amsq1,x1,ex1 = 0.;
+//calculate Bk(j1) for gaussian W(m1) or non zero Sigma
+
+	A[0] = j1; 
+	A[1] = j1; 
+	
+	double sfact = sqrt(2*j1 +1); 
+	double cgg = 0.; 
+	double tTerm = 0.; 
+	double II = 0.;
+	for(int i = 2; i<= 4; i = i + 2){
+		
+		A[2] = i; 
+		Bk11 = 0.;
+		Bk12 = 0.; 
+		for(int m = 0; m <= j14; m = m +2){
+			am11 = 0.5*(m - j12);
+			amsq1 = pow(am11,2);
+			x1 = -(amsq/(2*sigsq));
+			II = j1 - am11; 
+			A[3] = am11; 
+			A[4] = -am11; 
+			cgg = CG2(A);
+			ex1 = exp(x1);
+			tTerm = cn1*ex1*pow(-1,II) * sfact*cgg;
+			if(i == 2) Bk11 = Bk11 * tTerm;
+			if(i == 4) Bk12 = Bk12 * tTerm; 
+		}
+	}
+	
+
+//need to add correction of statistical tensors for cascade feeding. 
+//
+//
+
+
+//Calculate Rk(LL,j1j2)
+//
+	double B[6];
+	
+	double c0,c1,c2,q0,q1,q2,i0,i1,i2,sf0,sf1,sf2, rk01,rk11,rk21, rk02, rk12,rk22 = 0.; 
+	
+	for(int k = 1; k < 4; k = k + 2){
+
+		A[0] = AL0;
+		A[1] = AL0;
+		A[2] = k; 
+		A[3] = 1.;
+		A[4] = -1.;
+		A[5] = 0.;
+		
+		B[0] = j1; 
+		B[1] = j1; 
+		B[2] = AL0; 
+		B[3] = AL0; 
+		B[4] = k; 
+		B[5] = j2; 
+	
+		c0 = CG2(A);
+		q0 = RACAH(B);
+		i0 = 1+j1 -j2-k; 
+
+		sf0 = sqrt((2*j1+2) * (2*AL0 +1) * (2*AL1 +1)); 
+
+		if(k == 2){ rk01 = (pow(-1,i0) * sf0 * c0 * q0);}	
+		if(k == 4){ rk02 = (pow(-1,i0) * sf0 * c0 * q0);}
+
+		A[1] = AL1; 
+		B[3] = AL1; 
+
+		c1 = CG2(A);
+		q1 = RACAH(B);
+		i1 = 1 +j1 - j2 + AL1 - AL0 - k;
+		sf1 = sqrt((2*j1 +1) * (2*AL0 +1) *(2*AL1 +1));
+		 
+		if(k == 2){ rk11 = (pow(-1,i1) * sf1 * c1 * q1);}
+		if(k == 4){ rk12 = (pow(-1,i1) * sf1 * c1 * q1);}
+
+		A[0] = AL1; 
+		B[2] = AL1; 
+		
+		c2 = CG2(A); 
+		q2 = RACAH(B); 
+		i2 = 1 + j1 - j2 - k;
+		sf2 = sqrt((2*j1 +1) * (2*AL1+1) *(2*AL1+1));
+
+		if(k == 2){ rk21 = (pow(-1,i2) * sf2 * c2 * q2);}
+		if(k == 4){ rk22 = (pow(-1,i2) * sf2 * c2 * q2);}
+		
+
+	}
+
+
+
+
+
+	
+	printf("Bk1 = %lf\n",Bk11);
+	printf("Bk2 = %lf\n",Bk12);
+
+	printf("rk01 = %lf\n",rk01); 
+	printf("rk11 = %lf\n",rk11); 
+	printf("rk21 = %lf\n",rk21); 
+	printf("rk02 = %lf\n",rk02); 
+	printf("rk12 = %lf\n",rk12); 
+	printf("rk22 = %lf\n",rk22); 
+	
+
+/*
+
+	double JJ[6] = {j1,j1,2,.5,-.5,0};
+	double racah = RACAH(JJ);	
+	double cg222 = CG2(JJ);
+
+	
+	
+//	double cg = CG(J,M,j1,m1,j2,m2);
+//	double W3J= ThreeJsym(j1,m1,j2,m2,j3,m3);
+//	double W6J = SixJsym(j1,j2,j3,j4,j5,j6);
+
+	double cg_a = CG_a(JJ);
+	double W3J_a= ThreeJsym_a(JJ);
+	double W6J_a = SixJsym_a(JJ);
+	printf(" CG convert : %lf\n",cg222);	
+	printf(" CG         : %lf\n", cg_a);
+	printf(" W3-J 	    : %lf\n", W3J_a);
+	printf(" W6-J       : %lf\n", W6J_a);	
+	printf(" Racah convert  : %lf\n", racah);
+
+*/	
+	// here the Chi-sqs from A2 and A4 needed to be added together. 
+
+	 
+	double delta_min = -3.14159 / 2.;
+	double delta_max = 3.14159 / 2.;
+	double step = 0.001; 
+
+ 	double delta = 0.;
+	double tan_delta =0.;
+	double A0E = residual[0];//NEED FrOM AF FIT
+	double A2E = residual[1];
+	double A4E = residual[2]; 
+	double A2T,a2T = 0.;
+	double A4T,a4T = 0.;
+	double a4E = residual[2]; //NEED FROM AD FIT
+	double rd2T = 0.;
+	double rd4T = 0.;
+	double X22,X24 = 0.;
+	double X2_total = 0.;
+
+	int points = (delta_max - delta_min) / step ; 
+	
+	vector<double> chisqr;
+	vector<double> tdelta;
+
+	for(int i = 0; i< points; i++){
+		tan_delta = i + step + delta_min;
+		delta = tan(tan_delta);
+		rd2T = (rk01 + 2*delta*rk11 + pow(delta,2)*rk21)/(1+pow(delta,2));
+		A2T = QD2*Bk11*rd2T;
+		a2T = A2T/A0E;
+
+		
+		rd4T = (rk02 + 2*delta*rk12 + pow(delta,2)*rk22)/(1+pow(delta,2));
+		A4T = QD4*Bk12*rd4T;
+		a4T = A4T/A0E;
+	
+		X22 = pow(abs(A2E- a2T),2)/(abs(a2T));
+		X24 = pow(abs(A4E -a4T),2)/(abs(a4T));
+		X2_total = X22 + X24;
+
+		chisqr.push_back(X2_total);
+		tdelta.push_back(tan_delta);
+
+	}
+
+
+
+
+
+
+	int param;
+
+	param = param_run(det_param_token, gamma_energy_token, ang_file_token, sigma_token, feeding_token, j1j2token);  	
+
+
+	int optnum = -1;
+	menu();
+
+	printf("Input option : ");
+	scanf("%d", &optnum);
+	if(optnum < 0){
+		optnum = -1;
+		printf("Negative numbers are not allowed!\nInput option : ");
+		scanf("%d", &optnum);
+	}
+
+	if(optnum == 0){
+		Readme();
+		optnum = -1;
+
+		printf("Input option : ");
+		scanf("%d", &optnum);
+	} 
+
+	//if all inputs are valid, plot distribution. 
+	if(param == 1 && optnum == 1){
+		optnum = -1;
+
+		//plotting values here
+		std::vector<double> x;
+		std::vector<double> y;
+		
+	        //x will be arctan of mixing ratio. 
+	        //y will be log(X^2); 
+	        
+				
+
+		for(int i = 0; i < 50; i++){
+			x.push_back( i );
+			y.push_back( pow(i, 0.5) );
+		}
+		gui.SetData(x,y);
+		gui.Init();
+		gui.Loop();
+		gui.Close();
+/*
+		printf("Input option : ");
+		scanf("%d", &optnum);
+		if(optnum < 0){		
+			printf("Input option : ");
+			scanf("%d", &optnum);
+		}*/	
+	}
+
+
+
+	return 1;
+
+
+}
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Cleb.cxx b/Cleb.cxx
new file mode 100644
index 0000000..4c3796c
--- /dev/null
+++ b/Cleb.cxx
@@ -0,0 +1,740 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//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) = %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));
+					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;
+
+			}
+// 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;
+					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){
+	//recall that j1,m1 + j2,m2 = J,M
+
+	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 A = 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*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){
+
+
+	double A[6] = {1,1,2,0,0,0};
+
+
+	double cg = CG2(A);
+	double cgr = CG(1,1,2,0,0,0);
+	printf("------\n");
+	printf("CG = %lf\n",cg);
+	printf("CGR = %lf\n",cgr);
+
+
+	return 0; 
+}
+
+
+
diff --git a/Cleb.h b/Cleb.h
new file mode 100644
index 0000000..4d4be1e
--- /dev/null
+++ b/Cleb.h
@@ -0,0 +1,653 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//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;
+
+					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.
+					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;
+
+					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;
+
+			}
+// 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;
+
+					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.
+					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;
+
+					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){
+	//recall that j1,m1 + j2,m2 = J,M
+
+	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 A = 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*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;
+}
+
+
diff --git a/Coeff.h b/Coeff.h
new file mode 100644
index 0000000..de6ca97
--- /dev/null
+++ b/Coeff.h
@@ -0,0 +1,71 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
+
+using namespace std;
+
+
+int param_run(int dt, int gt, int at, int st, int ft, int jt){
+	int param = 0; 
+	if( dt == 1 && gt == 1 && at == 1 && st == 1 && ft == 1 && jt == 1){
+		param =0; 
+	}else param = 1; 
+
+	return param;
+}
+
+
+
+int StrToInt(std::string const& s)
+{
+    std::istringstream iss(s);
+    int value;
+
+    if (!(iss >> value)) throw std::runtime_error("invalid int");
+
+    return value;
+}
+
+
+int factorial(int fact){
+	for(int i=1;i<=fact;i++){
+		fact=fact*i;
+	} 
+	return fact;
+}
+
+
+double FACTLOG(int num){
+
+	double faclog[170];
+	int RI = 0;
+
+	faclog[0] = 0.;
+	faclog[1] = 0.;
+	
+	for( int i = 3; i < 170; i++){
+		RI = i - 1; 
+		faclog[i] = log(RI) + faclog[i-1];
+
+	}
+	
+	double flog = faclog[num];
+
+
+	return flog;
+}
+
+
+
+
+
+
+
diff --git a/GUI.h b/GUI.h
new file mode 100644
index 0000000..dea6673
--- /dev/null
+++ b/GUI.h
@@ -0,0 +1,422 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//To compile : g++ AD.cxx -o {Input Executable Name} -lX11
+// VERSION : 1.0
+// Created with thanks to Daniel Folds Holt - University of Cambridge Cavendish Lab
+
+
+
+
+using namespace std;
+
+
+class HistoGUI{
+
+
+	public:
+
+		HistoGUI(){}
+
+		Display * disp;
+		Window    wind;
+		XEvent    evt;
+		int       screen;
+
+		XColor xcolour;
+		XColor xcolour_red;
+		XColor xcolour_veryred;
+		Colormap cmap;
+
+
+		std::vector<double> x;
+		std::vector<double> y;
+		double max_x;
+		double min_x;
+		double max_y;
+		double min_y;
+
+		double pos_x;
+		double pos_y;
+		double new_pos_x;
+		double new_pos_y;
+
+		unsigned int width;  
+		unsigned int height;  
+		double window_width;
+		double window_height;
+		double width_scale; 
+		double x_offset;
+		double height_scale;
+		double y_offset;
+
+		double old_xl, old_xh, old_yl, old_yh;
+		double old_mouse_x, old_mouse_y;
+
+		int Init();
+		int SetData(std::vector<double> a, std::vector<double>b);
+		int Loop();	
+		void Close(){ printf("Closing now!\n"); XCloseDisplay(disp); }
+		int DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win);
+		int DrawCrosshairs(int mouse_x, int mouse_y);
+		int Zoom(int mouse_x, int mouse_y);
+
+
+//		int Legendre_fit(std::vector<double> d1, std::vector<double> d2, double 
+
+};
+
+int HistoGUI::SetData(std::vector<double> a, std::vector<double>b){
+	for(int i=0; i < a.size(); i++){
+		x.push_back(a[i]);
+		y.push_back(b[i]);
+	}
+	return a.size();
+}
+
+int HistoGUI::Init(){
+
+	disp = XOpenDisplay(NULL);
+	if (disp == NULL) {
+		fprintf(stderr, "Cannot open display\n");
+		exit(1);
+	}
+
+	screen = DefaultScreen(disp);
+	wind   = XCreateSimpleWindow(disp, RootWindow(disp, screen), 10, 10, 600, 400, 1,
+			BlackPixel(disp, screen), WhitePixel(disp, screen));
+	XGrabPointer(disp, wind, False, ButtonPressMask | ButtonReleaseMask | PointerMotionMask, GrabModeAsync,
+			GrabModeAsync, None, None, CurrentTime);
+	XSelectInput(disp, wind, ExposureMask | KeyPressMask | ButtonPressMask | ButtonReleaseMask | PointerMotionMask);
+	XMapWindow(disp, wind);
+
+	// colours
+	cmap = DefaultColormap(disp, screen);
+	xcolour.red = 32000; 
+	xcolour.green = 32000; 
+	xcolour.blue = 42000;
+	xcolour.flags = DoRed | DoGreen | DoBlue;
+	XAllocColor(disp, cmap, &xcolour);
+
+	xcolour_red.red   = 42000; 
+	xcolour_red.green = 32000; 
+	xcolour_red.blue  = 32000; 
+	xcolour_red.flags = DoRed | DoGreen | DoBlue;
+	XAllocColor(disp, cmap, &xcolour_red);
+
+	xcolour_veryred.red   = 65000; 
+	xcolour_veryred.green = 32000; 
+	xcolour_veryred.blue  = 32000; 
+	xcolour_veryred.flags = DoRed | DoGreen | DoBlue;
+	XAllocColor(disp, cmap, &xcolour_veryred);
+
+
+
+	return 1;
+
+}
+
+int HistoGUI::DrawCrosshairs(int mouse_x, int mouse_y){
+
+	int j1, j2;
+	unsigned int j3, j4;  
+	Window root_return;
+
+	XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4);
+	XSetForeground(disp,  DefaultGC(disp, screen), xcolour_red.pixel);
+	XDrawLine(disp, wind, DefaultGC(disp, screen), mouse_x, 0.0, mouse_x, height);
+	XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, mouse_y, width,  mouse_y);
+
+	pos_y = mouse_y * height_scale - y_offset;
+	pos_x = mouse_x * width_scale  - x_offset;
+
+	char coord[32];
+	sprintf(coord, "(%.2f, %.2f)", pos_x, pos_y);
+	XDrawString(disp, wind, DefaultGC(disp, screen), mouse_x + 5, mouse_y + 12, coord, strlen(coord));
+
+	return 1;
+}
+
+int HistoGUI::Zoom(int mouse_x, int mouse_y){
+
+	new_pos_x = mouse_x * width_scale - x_offset;
+	new_pos_y = mouse_y * height_scale  - y_offset;
+
+	pos_x = old_mouse_x * width_scale - x_offset;
+	pos_y = old_mouse_y * height_scale  - y_offset;
+
+
+	if(std::abs(mouse_x - old_mouse_x) < 20 or std::abs(mouse_y - old_mouse_y) < 20) return 1;
+
+
+	//printf("Zooming: %f - %f, %f - %f\n", mouse_x, old_mouse_x,mouse_y,old_mouse_y);
+
+	double x_low_win, y_low_win, x_hi_win, y_hi_win;
+	if(new_pos_x < pos_x){
+		x_low_win = new_pos_x;
+		x_hi_win  = pos_x;
+		//printf("new low; [%f, %f]\n",x_low_win, x_hi_win);
+	} else {
+		x_low_win = pos_x;
+		x_hi_win  = new_pos_x;
+		//printf("new high; [%f, %f]\n",x_low_win, x_hi_win);
+	}
+	if(new_pos_y > pos_y){
+		y_low_win = new_pos_y;
+		y_hi_win  = pos_y;
+	} else {
+		y_low_win = pos_y;
+		y_hi_win  = new_pos_y;
+	}
+
+
+	//printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win);
+	XClearWindow(disp, wind);
+	DrawData(x_low_win, y_low_win, x_hi_win, y_hi_win);
+	//double x_low_win, double y_low_win, double x_hi_win, double y_hi_win
+
+	return 1;
+
+}
+
+
+int HistoGUI::DrawData(double x_low_win, double y_low_win, double x_hi_win, double y_hi_win){
+
+
+	int j1, j2;
+	unsigned int j3, j4;  
+	Window root_return;
+
+	XGetGeometry(disp, wind, &root_return, &j1, &j2, &width, &height, &j3, &j4);
+
+	//printf("Width = %u, Height = %u, x = %i, y = %i\n", width, height, j1, j2);
+	//printf("[(%f,%f) (%f,%f)]\n",x_low_win, y_low_win, x_hi_win, y_hi_win);
+
+	double x_step;// = (width  * 0.8) / x.size();
+	double y_step;// = (height * 0.8) / x.size();
+
+
+	if(x_low_win == -1 and y_low_win == -1 and x_hi_win == -1 and y_hi_win == -1){	
+		// Audomatically decide data postition
+		max_x = x[0];
+		max_y = y[0];
+		min_x = x[0];
+		min_y = y[0];
+
+		for(int i=0; i<x.size(); i++){
+			if(x[i] > max_x) max_x = x[i];
+			if(x[i] < min_x) min_x = x[i];
+			if(y[i] > max_y) max_y = y[i];
+			if(y[i] < min_y) min_y = y[i];
+		}
+		//printf(" max_x = %f\n", max_x );
+		//printf(" max_y = %f\n", max_y );
+		//printf(" min_x = %f\n", min_x );
+		//printf(" min_y = %f\n", min_y );
+
+		width_scale = (max_x - min_x) / (0.8 * width);
+		x_offset = 0.5 * ((max_x - min_x) / 0.8) - 0.5 * (min_x + max_x);
+
+		height_scale = -1. * (max_y - min_y) / (0.8 * height);
+		y_offset = -0.5 * ((max_y - min_y) / 0.8) + -0.5 * (min_y + max_y);
+
+		//printf("width scale  = %f\n", width_scale);
+		//printf("x_offset     = %f\n", x_offset);
+		//printf("height scale = %f\n", height_scale);
+		//printf("y_offset     = %f\n", y_offset);
+		double axis_x  = (0. + x_offset) / width_scale;
+		double axis_y  = (0. + y_offset) / height_scale;
+
+		XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel);
+		XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height);
+		XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y);
+
+		char axis_val[4];
+		int w_step = width / 10;
+		for(int i=0; i < (int) width; i += w_step){
+			double x_val   = i * width_scale - x_offset;
+			sprintf(axis_val, "%.1f", x_val);
+			XDrawString(disp, wind, DefaultGC(disp, screen), i, axis_y + 10, axis_val, strlen(axis_val));
+		}
+
+		int h_step = height / 10;
+		for(int i=0; i < (int) height; i += h_step){
+			double y_val = i * height_scale - y_offset;
+			sprintf(axis_val, "%.1f", y_val);
+			XDrawString(disp, wind, DefaultGC(disp, screen), axis_x + 10, i, axis_val, strlen(axis_val));
+		}
+
+
+		XSetForeground(disp, DefaultGC(disp,screen), 0);	
+		double x_wid ;
+		double y_wid ;
+		double x_wid2;
+		double y_wid2;
+
+		for(int i=0; i < x.size() - 2; i++){
+			x_wid  = (x[i] + x_offset) / width_scale;
+			y_wid  = (y[i] + y_offset) / height_scale;
+			x_wid2 = (x[i + 1] + x_offset) / width_scale;
+			y_wid2 = (y[i + 1] + y_offset) / height_scale;
+			//printf("(%f, %f), (%f,%f)\n", x_wid,y_wid,x_wid2,y_wid2);
+			XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2);
+			XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4);
+		}
+		XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4);
+
+	} else {
+		//double x_low_win, double y_low_win, double x_hi_win, double y_hi_win
+
+		width_scale = (x_hi_win - x_low_win) / width;
+		x_offset = -1.0 * x_low_win;
+
+		height_scale = (y_hi_win - y_low_win) / height;
+		y_offset = -1.0 * y_low_win;
+
+		//printf("width scale  = %f\n", width_scale);
+		//printf("x_offset     = %f\n", x_offset);
+		//printf("height scale = %f\n", height_scale);
+		//printf("y_offset     = %f\n", y_offset);
+		double axis_x  = (0. + x_offset) / width_scale;
+		double axis_y  = (0. + y_offset) / height_scale;
+
+		XSetForeground(disp, DefaultGC(disp,screen), xcolour.pixel);
+		XDrawLine(disp, wind, DefaultGC(disp, screen), axis_x, 0.0, axis_x, height);
+		XDrawLine(disp, wind, DefaultGC(disp, screen), 0.0, axis_y, width, axis_y);
+
+		char axis_val[4];
+		int w_step = width / 10;
+		for(int i=0; i < (int) width; i += w_step){
+			double x_val   = i * width_scale - x_offset;
+			sprintf(axis_val, "%.1f", x_val);
+			XDrawString(disp, wind, DefaultGC(disp, screen), i, axis_y + 10, axis_val, strlen(axis_val));
+		}
+
+		int h_step = height / 10;
+		for(int i=0; i < (int) height; i += h_step){
+			double y_val = i * height_scale - y_offset;
+			sprintf(axis_val, "%.1f", y_val);
+			XDrawString(disp, wind, DefaultGC(disp, screen), axis_x + 10, i, axis_val, strlen(axis_val));
+		}
+
+
+		XSetForeground(disp, DefaultGC(disp,screen), 0);	
+		double x_wid ;
+		double y_wid ;
+		double x_wid2;
+		double y_wid2;
+
+		for(int i=0; i < x.size() - 2; i++){
+			x_wid  = (x[i] + x_offset) / width_scale;
+			y_wid  = (y[i] + y_offset) / height_scale;
+			x_wid2 = (x[i + 1] + x_offset) / width_scale;
+			y_wid2 = (y[i + 1] + y_offset) / height_scale;
+			//	printf("(%f, %f), (%f,%f)\n", x_wid,y_wid,x_wid2,y_wid2);
+			XDrawLine(disp, wind, DefaultGC(disp, screen), x_wid, y_wid, x_wid2, y_wid2);
+			XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid -2, y_wid -2, 4, 4);
+		}
+		XFillRectangle(disp, wind, DefaultGC(disp, screen), x_wid2 -2, y_wid2 -2, 4, 4);
+
+
+	}
+
+	old_xl = x_low_win;
+	old_yl = y_low_win;
+	old_xh = x_hi_win;
+	old_yh = y_hi_win;
+
+	return 1;
+}
+
+
+int HistoGUI::Loop(){
+
+	bool MousePressed  = false;
+	bool MousePressed2 = false;
+
+	while (1) {
+		XNextEvent(disp, &evt);
+		if (evt.type == Expose) {
+			//		XFillRectangle(disp, wind, DefaultGC(disp, screen), 20, 20, 10, 10);
+			//			XDrawString   (disp, wind, DefaultGC(disp, screen), 10, 50, msg, strlen(msg));
+			DrawData(-1,-1,-1,-1);
+		} else if (evt.type == ButtonPress){
+			/* store the mouse button coordinates in 'int' variables. */
+			/* also store the ID of the window on which the mouse was */
+			/* pressed.                                               */
+			int mouse_x = evt.xbutton.x;
+			int mouse_y = evt.xbutton.y;
+			old_mouse_x = mouse_x;
+			old_mouse_y = mouse_y;
+
+			/* check which mouse button was pressed, and act accordingly. */
+			if(evt.xbutton.button == Button1){
+				/* draw a pixel at the mouse position. */
+				XClearWindow(disp, wind);
+				DrawCrosshairs(mouse_x, mouse_y);
+				DrawData(old_xl, old_yl, old_xh, old_yh);
+				MousePressed = true;
+				printf("(%f, %f)\n", pos_x, pos_y);
+			} if(evt.xbutton.button == Button3){
+				MousePressed2 = true;
+			} if(evt.xbutton.button == Button2){
+				XClearWindow(disp, wind);
+				DrawData(old_xl, old_yl, old_xh, old_yh);
+			}
+		} else if (evt.type == ButtonRelease){
+			/* store the mouse button coordinates in 'int' variables. */
+			/* also store the ID of the window on which the mouse was */
+			/* pressed.                                               */
+			int mouse_x = evt.xbutton.x;
+			int mouse_y = evt.xbutton.y;
+			//printf("Released button: %i, %i\n", mouse_x, mouse_y);
+			//printf("Pressed button: %i, %i\n", mouse_x, mouse_y);
+			/* check which mouse button was pressed, and act accordingly. */
+			if(evt.xbutton.button == Button1){
+				/* draw a pixel at the mouse position. */
+				Zoom(mouse_x, mouse_y);
+				//DrawData();
+				MousePressed = false;
+			} if(evt.xbutton.button == Button3){
+				MousePressed2 = false;
+			}
+		} else if (evt.type == MotionNotify and MousePressed){
+
+			int mouse_x = evt.xmotion.x;
+			int mouse_y = evt.xmotion.y;
+			XClearWindow(disp, wind);
+
+			DrawCrosshairs(old_mouse_x, old_mouse_y);
+			DrawCrosshairs(mouse_x, mouse_y);
+			DrawData(old_xl, old_yl, old_xh, old_yh);
+
+		} else if (evt.type == MotionNotify and MousePressed2){
+
+			int mouse_x = evt.xmotion.x;
+			int mouse_y = evt.xmotion.y;
+
+			XSetForeground(disp, DefaultGC(disp,screen), xcolour_veryred.pixel);
+			XFillRectangle(disp, wind, DefaultGC(disp, screen), mouse_x -2, mouse_y -2, 4, 4);
+
+		} else if (evt.type == KeyPress){
+
+			printf("Key pressed = %x\n", evt.xkey.keycode);
+			if(evt.xkey.keycode == 0x41 or evt.xkey.keycode == 0x39){
+				XClearWindow(disp, wind);
+				DrawData(-1,-1,-1,-1);
+
+			} else {
+				break;	
+			}
+		}
+	}
+
+	return 1;
+}
diff --git a/QDK.h b/QDK.h
new file mode 100644
index 0000000..b915169
--- /dev/null
+++ b/QDK.h
@@ -0,0 +1,206 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+
+using namespace std;
+
+
+int QKN(double Energy, double radius, double distance, double thickness){
+	
+	
+	double Qkn = 0.; 
+
+	double E_mev = Energy/1000;
+	
+	double E_log = log(E_mev);
+
+	double EL1 = E_log;
+	double EL2 = pow(E_log,2);
+	double EL3 = EL1*EL2;
+	double EL4 = pow(EL2,2);
+	double EL5 = EL4*EL1;
+
+	double TT = -1.1907 -0.5372*EL1 - 0.0438*EL2 + 0.0218*EL3 + 0.0765*EL4 + 0.0095*EL5;
+
+	double Tau = exp(TT);
+
+	//calulating attenuation angles
+	
+	
+	double Z1 = radius / (distance + thickness);
+	double Z2 = radius / distance; 
+
+	double alpha = atan(Z1);
+	double gamma = atan(Z2);
+	double beta;
+
+
+	double BL = 0.;
+	double BU = alpha;
+	double A = 0.;
+	double delx1  = (BU-BL)/1000;
+	
+	double sum1,sum2,sum3 = 0.;
+	double sum4,sum5,sum6 = 0.;
+
+	double cosb,sinb,secb,c2,c4,fac1,ex1 = 0.;
+ 	double term1,term2,term3 = 0.;
+	double term4,term5,term6 = 0.;
+	int J=0;
+	
+	int loop_length = 500;
+
+	for(int i = 0; i<=loop_length; i++){
+		
+		if(i > 0 and i < loop_length){
+			J = i%2;	
+			//printf("\t\ti = %d\nJ=%d\n",i,J);
+			if(J==0){A=2.;
+			}else {A=4.;}
+			beta = BL+i+delx1;
+		}else{A=.1;beta = BL+i+delx1;}
+	
+		cosb = cos(beta);
+		sinb = sin(beta);
+		secb = 1.0/cosb;
+		c2 = pow(cosb,2);		
+		c4 = pow(cosb,4);
+		fac1 = -1 *Tau *thickness *secb;
+		ex1 = exp(fac1);
+
+		term1 = 0.5*(3*c2-1)*(1-ex1)*sinb*A*delx1;
+		term2 = 0.125*A*(35*c4-30*c2+3)*(1-ex1)*sinb*delx1;
+		term3 = A*(1-ex1)*sinb*delx1;
+
+		sum1 = sum1 +term1;
+		sum2 = sum2 +term2;
+		sum3 = sum3 +term3;
+
+	}
+		
+
+	double ans1 = sum1/3;
+	double ans2 = sum2/3;
+	double ans3 = sum3/3;	
+
+	double J2,B,beta2,cosb2,sinb2,secb2,cscb2,c22,c44,fac2,ex2 = 0.;
+	double LB=alpha;
+	double UB=gamma;
+	double delx2 = (UB-LB)/1000;
+
+	for(int i = 0; i<=loop_length; i++){
+		
+		if(i > 0 and i < loop_length){
+			J2 = i%2;	
+			//printf("\t\ti = %d\nJ=%d\n",i,J2);
+			if(J2==0){B=2.;
+			}else {B=4.;}
+			beta2 = LB+i+delx2;
+		}else{B=.1;beta2 = LB+i+delx2;}
+		
+		cosb2 = cos(beta2);
+		sinb2 = sin(beta2);
+		secb2 = 1.0/cosb2;
+		cscb2 = 1.0/sinb2;
+		c22 = pow(cosb2,2);		
+		c44 = pow(cosb2,4);
+		fac2 = -1 *Tau *(radius*cscb2 -distance*secb2);
+		ex2 = exp(fac2);
+/*
+		printf("cosb2---%lf\n",cosb2);
+		printf("sinb2---%lf\n",sinb2);
+		printf("secb2---%lf\n",secb2);
+		printf("cscb2---%lf\n",cscb2);
+		printf("c22---%lf\n",c22);
+		printf("c44---%lf\n",c44);
+		printf("fac2---%lf\n",fac2);
+		printf("ex2---%lf\n",ex2);
+		printf("-------------------------------------\n");
+*/
+		term4 = 0.5*(3*c22-1)*(1-ex2)*sinb2*B*delx2;
+		term5 = 0.125*B*(35*c44-30*c22+3)*(1-ex2)*sinb2*delx2;
+		term6 = B*(1-ex2)*sinb2*delx2;
+
+		sum4 = sum4 +term4;
+		sum5 = sum5 +term5;
+		sum6 = sum6 +term6;
+
+	}
+	double	ans4=sum4/3;
+	double	ans5=sum5/3;
+	double	ans6=sum6/3;
+
+/*
+	printf("ans1:%lf\n",ans1);
+	printf("ans2:%lf\n",ans2);
+	printf("ans3:%lf\n",ans3);
+	printf("ans4:%lf\n",ans4);
+	printf("ans5:%lf\n",ans5);
+	printf("ans6:%lf\n",ans6);
+*/	
+	double QD2  = (ans1+ans4)/(ans3+ans6);
+	double QD4  = (ans2+ans5)/(ans3+ans6);
+	printf("--------------\n");
+	printf("  QD2 = %lf\n",QD2);
+	printf("  QD4 = %lf\n",QD4);
+	printf("--------------\n");
+
+//Now output a file that contains R, D , T , gamma energy, attentuation coeff, q2 and q4
+
+	ofstream fileo;
+	fileo.open ("ad.txt");
+	fileo << "Radius = " << radius <<" [cm]\n";
+	fileo << "Distance = " << distance <<" [cm]\n";
+	fileo << "Thickness = " << thickness <<" [cm]\n";
+	fileo << "Atten.C = " << Tau <<" [cm^-1]\n";
+	fileo << "Gamma_E = " << Energy <<" [KeV]\n";
+	fileo << "QD2 = " << QD2 << "\n";
+	fileo << "QD4 = " << QD4 << "\n";
+	fileo.close();
+
+	return 1; 
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Racah.h b/Racah.h
new file mode 100644
index 0000000..2f7c85d
--- /dev/null
+++ b/Racah.h
@@ -0,0 +1,303 @@
+#include <X11/Xlib.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <vector>
+#include <cmath> 
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+//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;
+	
+	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;
+		}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{
+	
+	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;
+ 
+	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;
+}
+	
+	
+
diff --git a/ad.txt b/ad.txt
new file mode 100644
index 0000000..de93b13
--- /dev/null
+++ b/ad.txt
@@ -0,0 +1,7 @@
+Radius = 3 [cm]
+Distance = 4 [cm]
+Thickness = 5 [cm]
+Atten.C = 0.294297 [cm^-1]
+Gamma_E = 1062 [KeV]
+QD2 = -0.499975
+QD4 = 0.374938