KrazyHorse

17-02-2004, 01:36:50

#include <stdio.h>

#include <math.h>

#include "nr.h"

#include "nrutil.h"

main()

{

float *x, *y, *X, *Y;

float a, s0, M, g, w0;

float **s, **w;

float *Dx, *Dy;

float **f, **F;

float *I;

int i;

int j;

int N = 11;

x = vector (0, N-1);

y = vector (0, N-1);

X = vector (0, N-1);

Y = vector (0, N-1);

w = matrix (0, N-1, 0, N-1);

s = matrix (0, N-1, 0, N-1);

Dx = vector (0, N-1);

Dy = vector (0, N-1);

f = matrix (0, N-1, 0, N-1);

F = matrix (0, N-1, 0, N-1);

I = vector (0, N-1);

a = 0.054;

s0 = 8328.4;

M = 91.187;

g = 2.494;

w0 = 41.54;

for(i = 0; i < N; i++)

x[i] = i / 2;

for (j = 0; j < N; j++)

y[j] = j;

for (i = 0; i < N; i++)

X[i] = a * pow((1-x[i]), (a-1));

for (j = 0; j < N; j++)

Y[j] = a * pow((1-y[j]), (a-1));

for(i = 0; i < N; i++)

{for(j = 0; j < N; j++)

s[i][j] = x[i] * y[j] * s0;

}

for(i = 0; i < N; i++)

Dx[i] = a * pow((1 - x[i]), (a-1));

for(j = 0; j < N; j++)

Dy[j] = a * pow((1 - y[j]), (a-1));

for(i = 0; i < N; i++)

{for(j = 0; j < N; j++)

w[i][j] = w0 * s[i][j] * g * g / ( ( s[i][j] - M * M ) * ( s[i][j] - M * M ) + s[i][j] * s[i][j] * g * g / ( M * M));

}

for(i = 0; i < N; i++)

{for(j = 0; j < N; j++)

f[i][j] = Dx[i] * Dy[j] * w[i][j];

}

/* for(i = 0; i < N; i++) */

/* {for(j = 0; j < N; j++) */

/* f[i][j] = Dx[i] * Dy[j]; */

/* } */

F[0][j] = 0.0;

for(i = 2; i < N; i = i + 2)

{for(j = 0; j < N; j++)

F[i][j] = F[i-2][j] + (f[i-2][j] / 3 + 4 * f[i-1][j] / 3 + f[i][j] / 3) / (N-1);

}

I[0] = 0.0;

for(j = 2; j < N; j = j + 2)

I[j] = I[j-2] + (F[N-1][j-2] / 3 + 4 * F[N-1][j-1] / 3 + F[N-1][j] / 3) / (N-1);

/* for(i = 0; i < N; i= i + 2) */

printf("%18.15f\n", x[5]);

}

Now, most of this is irrelevant right now (it's supposed to be 2-d integration of some shit function via simpson's rule, but no matter)

I was getting a 0 for my integration, so I decided to do some checking.

The relevant part is

float *x, *y, *X, *Y;

int i;

int j;

int N = 11;

x = vector (0, N-1);

y = vector (0, N-1);

for(i = 0; i < N; i++)

x[i] = i / 2;

printf("%18.15f\n", x[5]);

I tell it to print 5/2 and it gives me 2.0000000

I defined x[i] as float; what's the problem?

#include <math.h>

#include "nr.h"

#include "nrutil.h"

main()

{

float *x, *y, *X, *Y;

float a, s0, M, g, w0;

float **s, **w;

float *Dx, *Dy;

float **f, **F;

float *I;

int i;

int j;

int N = 11;

x = vector (0, N-1);

y = vector (0, N-1);

X = vector (0, N-1);

Y = vector (0, N-1);

w = matrix (0, N-1, 0, N-1);

s = matrix (0, N-1, 0, N-1);

Dx = vector (0, N-1);

Dy = vector (0, N-1);

f = matrix (0, N-1, 0, N-1);

F = matrix (0, N-1, 0, N-1);

I = vector (0, N-1);

a = 0.054;

s0 = 8328.4;

M = 91.187;

g = 2.494;

w0 = 41.54;

for(i = 0; i < N; i++)

x[i] = i / 2;

for (j = 0; j < N; j++)

y[j] = j;

for (i = 0; i < N; i++)

X[i] = a * pow((1-x[i]), (a-1));

for (j = 0; j < N; j++)

Y[j] = a * pow((1-y[j]), (a-1));

for(i = 0; i < N; i++)

{for(j = 0; j < N; j++)

s[i][j] = x[i] * y[j] * s0;

}

for(i = 0; i < N; i++)

Dx[i] = a * pow((1 - x[i]), (a-1));

for(j = 0; j < N; j++)

Dy[j] = a * pow((1 - y[j]), (a-1));

for(i = 0; i < N; i++)

{for(j = 0; j < N; j++)

w[i][j] = w0 * s[i][j] * g * g / ( ( s[i][j] - M * M ) * ( s[i][j] - M * M ) + s[i][j] * s[i][j] * g * g / ( M * M));

}

for(i = 0; i < N; i++)

{for(j = 0; j < N; j++)

f[i][j] = Dx[i] * Dy[j] * w[i][j];

}

/* for(i = 0; i < N; i++) */

/* {for(j = 0; j < N; j++) */

/* f[i][j] = Dx[i] * Dy[j]; */

/* } */

F[0][j] = 0.0;

for(i = 2; i < N; i = i + 2)

{for(j = 0; j < N; j++)

F[i][j] = F[i-2][j] + (f[i-2][j] / 3 + 4 * f[i-1][j] / 3 + f[i][j] / 3) / (N-1);

}

I[0] = 0.0;

for(j = 2; j < N; j = j + 2)

I[j] = I[j-2] + (F[N-1][j-2] / 3 + 4 * F[N-1][j-1] / 3 + F[N-1][j] / 3) / (N-1);

/* for(i = 0; i < N; i= i + 2) */

printf("%18.15f\n", x[5]);

}

Now, most of this is irrelevant right now (it's supposed to be 2-d integration of some shit function via simpson's rule, but no matter)

I was getting a 0 for my integration, so I decided to do some checking.

The relevant part is

float *x, *y, *X, *Y;

int i;

int j;

int N = 11;

x = vector (0, N-1);

y = vector (0, N-1);

for(i = 0; i < N; i++)

x[i] = i / 2;

printf("%18.15f\n", x[5]);

I tell it to print 5/2 and it gives me 2.0000000

I defined x[i] as float; what's the problem?