Index: src/sas/models/c_extension/libigor/2Y_PairCorrelation.c
===================================================================
--- src/sas/models/c_extension/libigor/2Y_PairCorrelation.c	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
+++ src/sas/models/c_extension/libigor/2Y_PairCorrelation.c	(revision 794922220f2fb3b23dd08c648aabf87b4850596c)
@@ -0,0 +1,226 @@
+/*
+ *  PairCorrelation.c
+ *  twoyukawa
+ *
+ *  Created by Marcus Hennig on 5/9/10.
+ *  Copyright 2010 __MyCompanyName__. All rights reserved.
+ *
+ */
+
+#include "2Y_PairCorrelation.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+//#include <gsl/gsl_errno.h> 
+//#include <gsl/gsl_fft_real.h>
+
+/*
+=================================================================================================== 
+
+ Source: J.B.Hayter: A Program for the fast bi-directional transforms
+ between g(r) and S(Q), ILL internal scientific report, October 1979
+ 
+ The transformation between structure factor and pair correlation 
+ function is given by 
+ 
+ g(x) = 1 + 1 / (12*pi*phi*x) * int( [S(q)-q]*q*sin(q*x), { 0, inf } )
+ 
+ where phi is the volume frcation, x and q are dimensionless variables, 
+ scaled by the radius a of the particles:
+ 
+ r = x * a;
+ Q = q / a
+ 
+ Discretizing the integral leads to
+ 
+ x[k] = 2*pi*k / (N * dq)
+ g(x[k]) = 1 + N*dq^3 / (24*pi^2*phi*k) * Im{ sum(S[n]*exp(2*pi*i*n*k/N),{n,0,N-1}) }
+ 
+ where S[n] = n*(S(q[n])-1) with q[n]=n * dq
+
+=================================================================================================== 
+*/
+
+/*
+int PairCorrelation_GSL( double phi, double dq, double* Sq, double* dr, double* gr, int N )
+{
+	double* data = malloc( sizeof(double) * N );
+	int n;
+	for ( n = 0; n < N; n++ )
+		data[n] = n * ( Sq[n] - 1 );
+	
+	// data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 })
+	int stride = 1;
+	int error  = gsl_fft_real_radix2_transform( data, stride, N );
+	
+	// if no errors detected 
+	if ( error == GSL_SUCCESS ) 
+	{
+		double alpha = N * pow( dq, 3 ) / ( 24 * M_PI * M_PI * phi );
+	
+	
+		*dr = 2 * M_PI / ( N * dq );  
+		int k;
+		double real, imag;
+		for ( k = 0; k < N; k++ )
+		{
+			// the solutions of the transform is stored in data, 
+			// consult GSL manual for more details
+			if ( k == 0 || k == N / 2)
+			{ 
+				real = data[k];
+				imag = 0;
+			}
+			else if ( k < N / 2 ) 
+			{
+				real = data[k];
+				imag = data[N-k];
+			}
+			else if ( k > N / 2 )
+			{
+				real =  data[N-k];
+				imag = -data[k];	
+			}
+			
+			if ( k == 0 ) 
+				gr[k] = 0;
+			else
+				gr[k] = 1. + alpha / k * (-imag);
+		} 
+	} 
+	// if N is not a power of two
+	else if ( error == GSL_EDOM )
+	{
+		printf( "N is not a power of 2\n" );
+	}
+	else
+	{
+		printf( "Could not perform DFT (discrete fourier transform)\n" );
+	}
+	// release allocated memory
+	free( data );
+	
+	// return error value
+	return error;
+}
+ */
+
+
+// this uses numerical recipes for the FFT
+//
+int PairCorrelation( double phi, double dq, double* Sq, double* dr, double* gr, int N )
+{
+	double* data = malloc( sizeof(double) * N * 2);
+	int n,error,k;
+	double alpha,real,imag;
+	double Pi = 3.14159265358979323846264338327950288;   /* pi */
+
+	
+	for ( n = 0; n < N; n++ ) {
+		data[2*n] = n * ( Sq[n] - 1 );
+		data[2*n+1] = 0;
+	}
+	//	printf("start of new fft\n");
+	
+	// data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 })
+	//	int error  = gsl_fft_real_radix2_transform( data, stride, N );
+	error  = 1;
+	dfour1( data-1, N, 1 );		//N is the number of complex points
+	
+	//	printf("dfour1 is done\n");
+	
+	// if no errors detected 
+	if ( error == 1 ) 
+	{
+		alpha = N * pow( dq, 3 ) / ( 24 * Pi * Pi * phi );
+		
+		*dr = 2 * Pi / ( N * dq );  
+		for ( k = 0; k < N; k++ )
+		{
+			// the solutions of the transform is stored in data, 
+			// consult GSL manual for more details
+			if ( 2*k == 0 || 2*k == 2*N / 2)
+			{ 
+				real = data[2*k];
+				imag = 0;
+			}
+			else if ( 2*k < 2*N / 2 ) 
+			{
+				real = data[2*k];
+				imag = data[2*k+1];
+			}
+			else if ( 2*k > 2*N / 2 )
+			{
+				real =  data[2*k];
+				imag = -data[2*k+1];	
+			}
+			
+			if ( k == 0 ) 
+				gr[k] = 0;
+			else
+//				gr[k] = 1. + alpha / k * (-imag);		//if using GSL
+				gr[k] = 1. + alpha / k * (imag);		//if using NR
+		} 
+	} 
+	
+	// release allocated memory
+	free( data );
+	
+//	printf(" done with FFT assignment -- Using Numerical Recipes, not GSL\n");
+	
+	// return error value
+	return error;
+}
+
+
+// isign == 1 means no scaling of output. isign == -1 multiplies output by nn
+//
+//
+#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
+
+void dfour1(double data[], unsigned long nn, int isign)
+{
+	unsigned long n,mmax,m,j,istep,i;
+	double wtemp,wr,wpr,wpi,wi,theta;
+	double tempr,tempi;
+	
+	n=nn << 1;
+	j=1;
+	for (i=1;i<n;i+=2) {
+		if (j > i) {
+			SWAP(data[j],data[i]);
+			SWAP(data[j+1],data[i+1]);
+		}
+		m=n >> 1;
+		while (m >= 2 && j > m) {
+			j -= m;
+			m >>= 1;
+		}
+		j += m;
+	}
+	mmax=2;
+	while (n > mmax) {
+		istep=mmax << 1;
+		theta=isign*(6.28318530717959/mmax);
+		wtemp=sin(0.5*theta);
+		wpr = -2.0*wtemp*wtemp;
+		wpi=sin(theta);
+		wr=1.0;
+		wi=0.0;
+		for (m=1;m<mmax;m+=2) {
+			for (i=m;i<=n;i+=istep) {
+				j=i+mmax;
+				tempr=wr*data[j]-wi*data[j+1];
+				tempi=wr*data[j+1]+wi*data[j];
+				data[j]=data[i]-tempr;
+				data[j+1]=data[i+1]-tempi;
+				data[i] += tempr;
+				data[i+1] += tempi;
+			}
+			wr=(wtemp=wr)*wpr-wi*wpi+wr;
+			wi=wi*wpr+wtemp*wpi+wi;
+		}
+		mmax=istep;
+	}
+}
+#undef SWAP
