Problem with bessj0 in C program


dziadu
10-30-2006, 01:51 AM
Hi all!
It's my first post on this forum.
(sorry for my weak english)


I have to write three simple programs which interpolate points to Regular Spherical Bessel Function: one in Fortran wih NR, two in C with NR and GSL.

Programs in Fortran and C with GSL works fine, but wrote in C with NR don't work correctly:

There is source code of my program (03_lalik_1_nr.c):
// filename: 03_lalik_1_nr.c
#include <stdio.h>
#include <stdlib.h>

#define N 10

int main() {
float *x, *y, yy, dy;
int i;

x = (float*) malloc((N+1)*sizeof(float));
y = (float*) malloc((N+1)*sizeof(float));
// Preparing data...
for (i=1; i <= N; i++) {
x[i] = (10.0/(float)N)*i;
y[i] = bessj0(x[i]);
printf("%d: %f - %f\n", i, x[i], y[i]);
}

for (i=1; i <= N; i++) {
// Calculate interpolate points and print results
polint(x, y, N, x[i], &yy, &dy);
printf("%f\t%f\t%f\n", x[i], y[i], yy);
}
}

and this is result of execute:
1: 1.000000 - 134479872.000000
2: 2.000000 - 134494240.000000
3: 3.000000 - 134492192.000000
4: 4.000000 - 134490144.000000
5: 5.000000 - 134488096.000000
6: 6.000000 - 134486048.000000
7: nan - 134484000.000000
8: nan - 134482016.000000
9: nan - 134479968.000000
10: nan - 134494304.000000
Segmentation fault

This is my makefile
CC = gcc
CFLAGS =
LDFLAGS = -lm
RM = rm -f
GP = gnuplot

# make program
all: 03_lalik_1_nr

%.o : %.c
$(CC) -c $< -o $@

03_lalik_1_nr : bessj0.o 03_lalik_1_nr.o polint.o nrutil.o
$(CC) $(LDFLAGS) $^ -o $@

# draw
plot: 03_lalik_1
./$< > 03_lalik_1.outp
$(GP) 03_lalik_1.plot

clean:
$(RM) *.o 03_lalik_1.outp 03_lalik_1 03_lalik_1.png

help:
@echo "Usage:"
@echo "make help - pomoc"
@echo "make - generuje program wyjsciowy"
@echo "make plot - generuje program, wylicza dane oraz rysuje wykres"
@echo "make clean - czysci wszystkie smieci"

# depends
polint.o : polint.c
bessj0.o : bessj0.c
nrutil.o : nrutil.c
03_lalik_1_nr.o : 03_lalik_1_nr.c

I'll be appreciate for help.

All the best,
Rafał Lalik

..:: edit by dziadu

Ok, i've done some changes in sources, and now it looks like that:
// filnam: 03_lalik_1_nr.c
#include <stdio.h>
#include <stdlib.h>

#define N 10

int main() {
float *x, *y, yy, dy;
int i;
yy =0;
dy = 0;

x = vector(1, N);
y = vector(1, N);
// Preparing data...
for (i=0; i < N; i++) {
x[i] = (10.0/(float)N)*i;
y[i] = bessj0((float)x[i]);
printf("%d: %f - %f\n", i, x[i], y[i]);

}

for (i=0; i < N; i++) {
// Calculate interpolate points
polint(x, y, N, x[i], &yy, &dy);
// Print values
printf("%f\t%f\t%f\n", x[i], y[i], yy);
}
return 0;
}


I have tried to compile it and run on two machines. This is result:
AMD Athlon(tm) 64 Processor 3000+
x86_64-pc-linux-gnu-gcc-4.1.1
0: 0.000000 - 1065353216.000000
1: 1.000000 - 1065353216.000000
2: 2.000000 - 1065353216.000000
3: 3.000000 - 1065353216.000000
4: 4.000000 - 1065353216.000000
5: 5.000000 - 1065353216.000000
6: 6.000000 - 1065353216.000000
7: 7.000000 - 1065353216.000000
8: 8.000000 - 1065353216.000000
9: 9.000000 - 1065353216.000000
0.000000 1065353216.000000 0.000000
1.000000 1065353216.000000 0.000000
2.000000 1065353216.000000 0.000000
3.000000 1065353216.000000 0.000000
4.000000 1065353216.000000 0.000000
5.000000 1065353216.000000 0.000000
6.000000 1065353216.000000 0.000000
7.000000 1065353216.000000 0.000000
8.000000 1065353216.000000 0.000000
9.000000 1065353216.000000 0.000000

and Intel(R) Xeon(TM) CPU 3.00GHz
i485-pc-linux-gnu-gcc-3.4.4
0: 0.000000 - 1065353216.000000
1: 1.000000 - 1065353216.000000
2: 2.000000 - 1065353216.000000
3: 3.000000 - 1065353216.000000
4: 4.000000 - 1065353216.000000
5: 5.000000 - 1065353216.000000
6: 6.000000 - 1065353216.000000
7: 7.000000 - -4194304.000000
8: 8.000000 - -4194304.000000
9: 9.000000 - -4194304.000000
Segmentation fault


Please, please, help me!"