[230f479] | 1 | /* |
---|
| 2 | * PairCorrelation.c |
---|
| 3 | * twoyukawa |
---|
| 4 | * |
---|
| 5 | * Created by Marcus Hennig on 5/9/10. |
---|
| 6 | * Copyright 2010 __MyCompanyName__. All rights reserved. |
---|
| 7 | * |
---|
| 8 | */ |
---|
| 9 | |
---|
| 10 | #include "2Y_PairCorrelation.h" |
---|
| 11 | #include <stdio.h> |
---|
| 12 | #include <stdlib.h> |
---|
| 13 | #include <math.h> |
---|
| 14 | //#include <gsl/gsl_errno.h> |
---|
| 15 | //#include <gsl/gsl_fft_real.h> |
---|
| 16 | |
---|
| 17 | /* |
---|
| 18 | =================================================================================================== |
---|
| 19 | |
---|
| 20 | Source: J.B.Hayter: A Program for the fast bi-directional transforms |
---|
| 21 | between g(r) and S(Q), ILL internal scientific report, October 1979 |
---|
| 22 | |
---|
| 23 | The transformation between structure factor and pair correlation |
---|
| 24 | function is given by |
---|
| 25 | |
---|
| 26 | g(x) = 1 + 1 / (12*pi*phi*x) * int( [S(q)-q]*q*sin(q*x), { 0, inf } ) |
---|
| 27 | |
---|
| 28 | where phi is the volume frcation, x and q are dimensionless variables, |
---|
| 29 | scaled by the radius a of the particles: |
---|
| 30 | |
---|
| 31 | r = x * a; |
---|
| 32 | Q = q / a |
---|
| 33 | |
---|
| 34 | Discretizing the integral leads to |
---|
| 35 | |
---|
| 36 | x[k] = 2*pi*k / (N * dq) |
---|
| 37 | 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}) } |
---|
| 38 | |
---|
| 39 | where S[n] = n*(S(q[n])-1) with q[n]=n * dq |
---|
| 40 | |
---|
| 41 | =================================================================================================== |
---|
| 42 | */ |
---|
| 43 | |
---|
| 44 | /* |
---|
| 45 | int PairCorrelation_GSL( double phi, double dq, double* Sq, double* dr, double* gr, int N ) |
---|
| 46 | { |
---|
| 47 | double* data = malloc( sizeof(double) * N ); |
---|
| 48 | int n; |
---|
| 49 | for ( n = 0; n < N; n++ ) |
---|
| 50 | data[n] = n * ( Sq[n] - 1 ); |
---|
| 51 | |
---|
| 52 | // data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 }) |
---|
| 53 | int stride = 1; |
---|
| 54 | int error = gsl_fft_real_radix2_transform( data, stride, N ); |
---|
| 55 | |
---|
| 56 | // if no errors detected |
---|
| 57 | if ( error == GSL_SUCCESS ) |
---|
| 58 | { |
---|
| 59 | double alpha = N * pow( dq, 3 ) / ( 24 * M_PI * M_PI * phi ); |
---|
| 60 | |
---|
| 61 | |
---|
| 62 | *dr = 2 * M_PI / ( N * dq ); |
---|
| 63 | int k; |
---|
| 64 | double real, imag; |
---|
| 65 | for ( k = 0; k < N; k++ ) |
---|
| 66 | { |
---|
| 67 | // the solutions of the transform is stored in data, |
---|
| 68 | // consult GSL manual for more details |
---|
| 69 | if ( k == 0 || k == N / 2) |
---|
| 70 | { |
---|
| 71 | real = data[k]; |
---|
| 72 | imag = 0; |
---|
| 73 | } |
---|
| 74 | else if ( k < N / 2 ) |
---|
| 75 | { |
---|
| 76 | real = data[k]; |
---|
| 77 | imag = data[N-k]; |
---|
| 78 | } |
---|
| 79 | else if ( k > N / 2 ) |
---|
| 80 | { |
---|
| 81 | real = data[N-k]; |
---|
| 82 | imag = -data[k]; |
---|
| 83 | } |
---|
| 84 | |
---|
| 85 | if ( k == 0 ) |
---|
| 86 | gr[k] = 0; |
---|
| 87 | else |
---|
| 88 | gr[k] = 1. + alpha / k * (-imag); |
---|
| 89 | } |
---|
| 90 | } |
---|
| 91 | // if N is not a power of two |
---|
| 92 | else if ( error == GSL_EDOM ) |
---|
| 93 | { |
---|
| 94 | printf( "N is not a power of 2\n" ); |
---|
| 95 | } |
---|
| 96 | else |
---|
| 97 | { |
---|
| 98 | printf( "Could not perform DFT (discrete fourier transform)\n" ); |
---|
| 99 | } |
---|
| 100 | // release allocated memory |
---|
| 101 | free( data ); |
---|
| 102 | |
---|
| 103 | // return error value |
---|
| 104 | return error; |
---|
| 105 | } |
---|
| 106 | */ |
---|
| 107 | |
---|
| 108 | |
---|
| 109 | // this uses numerical recipes for the FFT |
---|
| 110 | // |
---|
| 111 | int PairCorrelation( double phi, double dq, double* Sq, double* dr, double* gr, int N ) |
---|
| 112 | { |
---|
| 113 | double* data = malloc( sizeof(double) * N * 2); |
---|
| 114 | int n,error,k; |
---|
| 115 | double alpha,real,imag; |
---|
| 116 | double Pi = 3.14159265358979323846264338327950288; /* pi */ |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | for ( n = 0; n < N; n++ ) { |
---|
| 120 | data[2*n] = n * ( Sq[n] - 1 ); |
---|
| 121 | data[2*n+1] = 0; |
---|
| 122 | } |
---|
| 123 | // printf("start of new fft\n"); |
---|
| 124 | |
---|
| 125 | // data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 }) |
---|
| 126 | // int error = gsl_fft_real_radix2_transform( data, stride, N ); |
---|
| 127 | error = 1; |
---|
| 128 | dfour1( data-1, N, 1 ); //N is the number of complex points |
---|
| 129 | |
---|
| 130 | // printf("dfour1 is done\n"); |
---|
| 131 | |
---|
| 132 | // if no errors detected |
---|
| 133 | if ( error == 1 ) |
---|
| 134 | { |
---|
| 135 | alpha = N * pow( dq, 3 ) / ( 24 * Pi * Pi * phi ); |
---|
| 136 | |
---|
| 137 | *dr = 2 * Pi / ( N * dq ); |
---|
| 138 | for ( k = 0; k < N; k++ ) |
---|
| 139 | { |
---|
| 140 | // the solutions of the transform is stored in data, |
---|
| 141 | // consult GSL manual for more details |
---|
| 142 | if ( 2*k == 0 || 2*k == 2*N / 2) |
---|
| 143 | { |
---|
| 144 | real = data[2*k]; |
---|
| 145 | imag = 0; |
---|
| 146 | } |
---|
| 147 | else if ( 2*k < 2*N / 2 ) |
---|
| 148 | { |
---|
| 149 | real = data[2*k]; |
---|
| 150 | imag = data[2*k+1]; |
---|
| 151 | } |
---|
| 152 | else if ( 2*k > 2*N / 2 ) |
---|
| 153 | { |
---|
| 154 | real = data[2*k]; |
---|
| 155 | imag = -data[2*k+1]; |
---|
| 156 | } |
---|
| 157 | |
---|
| 158 | if ( k == 0 ) |
---|
| 159 | gr[k] = 0; |
---|
| 160 | else |
---|
| 161 | // gr[k] = 1. + alpha / k * (-imag); //if using GSL |
---|
| 162 | gr[k] = 1. + alpha / k * (imag); //if using NR |
---|
| 163 | } |
---|
| 164 | } |
---|
| 165 | |
---|
| 166 | // release allocated memory |
---|
| 167 | free( data ); |
---|
| 168 | |
---|
| 169 | // printf(" done with FFT assignment -- Using Numerical Recipes, not GSL\n"); |
---|
| 170 | |
---|
| 171 | // return error value |
---|
| 172 | return error; |
---|
| 173 | } |
---|
| 174 | |
---|
| 175 | |
---|
| 176 | // isign == 1 means no scaling of output. isign == -1 multiplies output by nn |
---|
| 177 | // |
---|
| 178 | // |
---|
| 179 | #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr |
---|
| 180 | |
---|
| 181 | void dfour1(double data[], unsigned long nn, int isign) |
---|
| 182 | { |
---|
| 183 | unsigned long n,mmax,m,j,istep,i; |
---|
| 184 | double wtemp,wr,wpr,wpi,wi,theta; |
---|
| 185 | double tempr,tempi; |
---|
| 186 | |
---|
| 187 | n=nn << 1; |
---|
| 188 | j=1; |
---|
| 189 | for (i=1;i<n;i+=2) { |
---|
| 190 | if (j > i) { |
---|
| 191 | SWAP(data[j],data[i]); |
---|
| 192 | SWAP(data[j+1],data[i+1]); |
---|
| 193 | } |
---|
| 194 | m=n >> 1; |
---|
| 195 | while (m >= 2 && j > m) { |
---|
| 196 | j -= m; |
---|
| 197 | m >>= 1; |
---|
| 198 | } |
---|
| 199 | j += m; |
---|
| 200 | } |
---|
| 201 | mmax=2; |
---|
| 202 | while (n > mmax) { |
---|
| 203 | istep=mmax << 1; |
---|
| 204 | theta=isign*(6.28318530717959/mmax); |
---|
| 205 | wtemp=sin(0.5*theta); |
---|
| 206 | wpr = -2.0*wtemp*wtemp; |
---|
| 207 | wpi=sin(theta); |
---|
| 208 | wr=1.0; |
---|
| 209 | wi=0.0; |
---|
| 210 | for (m=1;m<mmax;m+=2) { |
---|
| 211 | for (i=m;i<=n;i+=istep) { |
---|
| 212 | j=i+mmax; |
---|
| 213 | tempr=wr*data[j]-wi*data[j+1]; |
---|
| 214 | tempi=wr*data[j+1]+wi*data[j]; |
---|
| 215 | data[j]=data[i]-tempr; |
---|
| 216 | data[j+1]=data[i+1]-tempi; |
---|
| 217 | data[i] += tempr; |
---|
| 218 | data[i+1] += tempi; |
---|
| 219 | } |
---|
| 220 | wr=(wtemp=wr)*wpr-wi*wpi+wr; |
---|
| 221 | wi=wi*wpr+wtemp*wpi+wi; |
---|
| 222 | } |
---|
| 223 | mmax=istep; |
---|
| 224 | } |
---|
| 225 | } |
---|
| 226 | #undef SWAP |
---|