Hopefully this works....
Code:
/*================================================= ==================*/
/* Program 'test' */
/* */
/* This program computes the self-diffusion on a lattice square */
/*-------------------------------------------------------------------*/
/* Author: Richard van Beusichem */
/* Email: [email protected] */
/*================================================= ==================*/
/*===== Include files ===============================================*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <nrutil.h>
#define ROW 5
#define COL 5
#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define M 7
#define NSTACK 50
/*===== Function prototypes =========================================*/
double expon(double x); // Returns an exponential random variable
double Random(); // Returns an uniformly distributed random number
int rand_int(int a, int b); // Returns random integer from interval [a,b]
void sort2(unsigned long n, float arr[], float brr[]);
/* Define the number of molecules */
const int MOL= 5;
/*====== Main Program ===============================================*/
int main() {
int i,e,a,b,lattice[ROW][COL]; /* array that represents the lattice */
double residence[ROW][COL]; /* array that holds exponential residence times */
double res_mol[MOL][2];
float mol_nr[MOL];
float time[MOL];
double exp_rv; /* exponential random variable */
double tau; /* mean residence time */
int x,y; /* x,y coordinates on lattice */
int j;
/* Initialize array */
for (i=0; i<COL; i++)
{
for (e=0; e<ROW; e++)
lattice[i][e]=0;
}
/* Initialize array */
for (i=0; i<COL; i++)
{
for (e=0; e<ROW; e++)
residence[i][e]=0;
}
/* Initialize array */
for (i=0; i<MOL; i++)
{
for (e=0; e<2; e++)
res_mol[i][e]=0;
}
a=0;
b=ROW-1;
tau=10.0;
/* Place molecules randomly on the lattice */
for (i=1;i<=MOL; i++)
{
x=rand_int(a,b);
y=rand_int(a,b);
while (lattice[x][y]>0)
{
x=rand_int(a,b);
y=rand_int(a,b);
}
lattice[x][y]=i;
}
for (i=0;i<ROW;i++)
{ for(e=0;e<COL;e++)
printf("%3d",lattice[i][e]);
printf("\n");
}
/* Generate and output exponential random variables */
for (i=0; i<ROW; i++)
{
for (e=0; e<COL; e++)
{
exp_rv = expon(1.0 / tau);
residence[i][e]=exp_rv;
}
}
for (i=0; i<ROW; i++)
{
for (e=0; e<COL; e++)
{
for (j=0; j<MOL; j++)
{
if (lattice[i][e]==j+1)
{
res_mol[j][0]=j+1;
res_mol[j][1]=residence[i][e];
}
}
}
}
for (i=0; i<MOL; i++)
{
mol_nr[i]=res_mol[i][0];
}
for (i=0; i<MOL; i++)
{time[i]=res_mol[i][1];
}
sort2(MOL,mol_nr,time);
for (i=0;i<ROW;i++)
{ for(e=0;e<COL;e++)
printf("%f",residence[i][e]);
printf("\n");
}
for (i=0;i<MOL;i++)
{ for(e=0;e<2;e++)
printf("%f",res_mol[i][e]);
printf("\n");
}
for (i=0; i<MOL; i++)
{
printf("%f",mol_nr[i]);
}
system("PAUSE");
return 0;
}
//================================================== =========================
//= Function to generate Exponentially Distributed Random Variables =
//= - Input: Mean value of distribution =
//= - Output: Returns with exponentially distributed random variable =
//================================================== =========================
double expon(double x)
{
double exp_value; // Computed exponential value to be returned
// Compute exponential random variable using inversion method
exp_value = -x * log(Random());
return(exp_value);
}
//================================================== =========================
//= Function to generate a random integer on a specified interval =
//= - Input: Interval [a,b] =
//= - Output: Returns with random integer from the interval =
//================================================== =========================
int rand_int(int a, int b)
{
return rand()%(b-a+1) + a;
}
//================================================== =========================
//= Function to generate Uniformly Distributed Random Numbers =
//= =
//= - Output: Returns uniformly distributed random numbers =
//================================================== =========================
static unsigned long int uliSeed = 1;
void InitRandom( long int iMySeed ) {
uliSeed = iMySeed;
}
double Random() {
long int liA = 16807;
long int liM = 2147483647;
long int liQ = 127773;
long int liR = 2836;
long int liTempSeed;
liTempSeed = liA*( uliSeed%liQ ) - liR*( uliSeed/liQ );
if( liTempSeed >= 0 )
uliSeed = liTempSeed;
else
uliSeed = liTempSeed + liM;
return( ((double) uliSeed)/liM );
}
//================================================== =========================
//= Function to sort an array arr[1..n] into ascending order =
//= using Quicksort, while making the corresponding rearrangement =
//= of the array brr[1..n]. =
//= =
//= - Input: Arrays =
//= - Output: Arrays sorted based on exponential residence times =
//================================================== =========================
void sort2(unsigned long n, float arr[], float brr[])
#include <nrutil.h>
{
unsigned long i,ir=n,j,k,l=1,*istack;
int jstack=0;
float a,b,temp;
istack=lvector(1,NSTACK);
for (;;) { /* Insertion sort when subarray small enough. */
if (ir-l < M) {
for (j=l+1;j<=ir;j++) {
a=arr[j];
b=brr[j];
for (i=j-1;i>=l;i--) {
if (arr[i] <= a) break;
arr[i+1]=arr[i];
brr[i+1]=brr[i];
}
arr[i+1]=a;
brr[i+1]=b;
}
if (!jstack) {
free_lvector(istack,1,NSTACK);
return;
}
ir=istack[jstack]; /* Pop stack and begin a new round of partitioning.*/
l=istack[jstack-1];
jstack -= 2;
} else {
k=(l+ir) >> 1; /* Choose median of left, center and right elements
as partitioning element a. Also
rearrange so that a[l] = a[l+1] = a[ir]. */
SWAP(arr[k],arr[l+1])
SWAP(brr[k],brr[l+1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
SWAP(brr[l],brr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
SWAP(brr[l+1],brr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
SWAP(brr[l],brr[l+1])
}
i=l+1; /*Initialize pointers for partitioning.*/
j=ir;
a=arr[l+1]; /*Partitioning element.*/
b=brr[l+1];
for (;;) { /*Beginning of innermost loop.*/
do i++; while (arr[i] < a); /*Scan up to find element > a.*/
do j--; while (arr[j] > a); /*Scan down to find element < a.*/
if (j < i) break; /*Pointers crossed. Partitioning complete.*/
SWAP(arr[i],arr[j]) /*Exchange elements of both arrays.*/
SWAP(brr[i],brr[j])
} /*End of innermost loop.*/
arr[l+1]=arr[j]; /*Insert partitioning element in both arrays.*/
arr[j]=a;
brr[l+1]=brr[j];
brr[j]=b;
jstack += 2;
/*Push pointers to larger subarray on stack, process smaller subarray immediately.*/
if (jstack > NSTACK) nrerror("NSTACK too small in sort2.");
if (ir-i+1 >= j-l) {
istack[jstack]=ir;
istack[jstack-1]=i;
ir=j-1;
} else {
istack[jstack]=j-1;
istack[jstack-1]=l;
l=i;
}
}
}
}
Below used header nrutil.h
Code:
void nrerror(char error_text[]);
float *vector(long nl, long nh);
int *ivector(long nl,long nh);
unsigned long *lvector(long nl,long nh);
double *dvector(long nl, long nh);
double **dmatrix(long nrl, long nrh, long ncl, long nch);
void free_vector(float *v, long nl, long nh);
void free_ivector(int *v, long nl,long nh);
void free_lvector(unsigned long *v, long nl,long nh);
void free_dvector(double *v, long nl, long nh);
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch);