source: sasview/sansmodels/src/libigor/2Y_PairCorrelation.c @ b11e127

/*
 *  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