Hi, I'm working on a program that models the trajectory of a baseball with drag. The program compiles okay and links okay. However, during runtime I get a segmentation fault. I remember when I was learning how to program in C several years ago that segmentation faults occur a lot in forgetting the & in scanf statments and when arrays access memory out of bounds. With this in mind, i took a look over the program and I didn't use any scanfs and arrays. I am using pointers though in which I'm not to used to. So I'm not sure how to track down where the segmentation fault is coming from.
I don't have a debugger since I'm compiling in cygwin so I'm pretty much on my own to figure where the problem exists. I was wondering if there is something obvious that any of you might see. I've pasted the code below for reference. Any suggestions would be greatly appreciated.
Thanks,
windell747
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define Tmax 500 // seconds
#define PI M_PI
/* this structure type for 3-vectors */
typedef struct{
double x, y, z;
} vec3;
/* we use a global time variable so the functions can use time
dependent accelerations if necessary*/
double t,t0, theta0;
double dragc(double H);
int main(int argc,char *argv[])
{
printf("hello\n");
vec3 xt,vt,xtold,vtold; //establish variables as vec3 types.
double dt;
double vt0,v0;
vec3 vecFRK2xv(), f_x(), f_v(), vec3sum(); //establish these functions as vec3 datatypes.
FILE *fout;
/*
if(argc<2){
fprintf(stderr, "usage: projectile [theta0 (deg)][v0 (m/s)]\n");
exit(0);
}
*/
fout = fopen("proj1.dat","w");
theta0 = (PI/180.0)*atof(argv[1]); // initial angle, radians
v0 = atof(argv[2]); // initial velocity, m/s
printf("theta0: %f, v0: %f\n",theta0,v0);
vtold.x = v0*cos(theta0); // x and y components of velocity
vtold.y = v0*sin(theta0);
vtold.z = 0.0; // in plane motion only here
t0 = 0.0; // initial time
xtold.x = 0.0; // initial position at origin
xtold.y = 0.0;
xtold.z = 0.0;
dt = 0.01; // time step
// starting values
fprintf(fout,"%e\t%e\t%e\t%e\t%e\n",t0, xtold.x,xtold.y,vtold.x,vtold.y);
for(t=t0; t<Tmax; t+= dt){
xt = vec3sum( xtold , vecFRK2xv(0,f_x,t,xtold,vtold,dt) );
vt = vec3sum( vtold , vecFRK2xv(1,f_v,t,xtold,vtold,dt) );
xtold = xt;
vtold = vt;
if(xt.y<0.0){
fprintf(fout,"%e\t%e\t%e\t%e\t%e\n", t,xt.x,xt.y,vt.x,vt.y);
break;
}
fprintf(fout,"%e\t%e\t%e\t%e\t%e\n", t,xt.x,xt.y,vt.x,vt.y);
}
fclose(fout);
}
// this function is complete: gives the formal definition dx/dt=v
vec3 f_x(double t, vec3 x, vec3 v)
{
return(v);
}
// this function should return the force/mass=acceleration as a
// function of the position and velocity vectors, based on the force law
vec3 f_v(double t, vec3 x, vec3 v)
{
vec3 f;
double b,h,magv,m,g;
h=x.y;
m=0.145; //mass in kilograms.
g=9.81; //we are assuming here that the gravitational constant is constant.
magv=sqrt((x.x*x.x)+(x.y*x.y));
f.x=-dragc(h)*magv*v.x;
f.y=-(m*g)-(dragc(h)*magv*v.y);
return(f);
}
// this function should return the drag coeffcient as a function of altitude H
// accounting for the density of air using an exponential scale height
double dragc(double H)
{
double b,h0,Cd,rho0,A;
h0=8300.0; //scale height of atmosphere.
Cd=0.35;
rho0=1.25; //density of air at sea level. kg/m^3
A=42.7762/10000.0; //cross-sectional area of baseball.
b=(0.5)*Cd*A*rho0*exp(-H/h0);
return(b);
}
typedef struct{
double x, y, z;
} vec3;
vec3 vecFRK2xv(int ytype, vec3 (*f_x)(double, vec3, vec3),
vec3 (*f_v)(double, vec3, vec3),
double t,vec3 xold,vec3 vold,double dt)
/*xold, vold, previous position and velocity
dt = step size
vec3 (*f_xv)(vec3, vec3) ==> pointer to a function of
two vec3s that returns a vec3. Either dv/dt or dx/dt.
*/
{
vec3 vec3sum(),scalar_vec3sum(),scalar_vec3mult();
vec3 k1x, k1v, k2;
k1x = scalar_vec3mult( dt, (*f_x)(t, xold, vold));
k1v = scalar_vec3mult( dt, (*f_v)(t, xold, vold));
if(ytype==0){
k2 = scalar_vec3mult( dt,
(*f_x)( t+dt/2.0, vec3sum(xold,scalar_vec3mult(0.5,k1x)),
vec3sum(vold,scalar_vec3mult(0.5,k1v))));
}else{ /* ytype = 1 */
k2 = scalar_vec3mult( dt,
(*f_v)( t+dt/2.0, vec3sum(xold,scalar_vec3mult(0.5,k1x)),
vec3sum(vold,scalar_vec3mult(0.5,k1v))));
}
return(k2); // add this value to the previous v or x
}
/* these are general vector routines needed by above */
vec3 vec3sum( vec3 X, vec3 Y)
{
vec3 tmp;
tmp.x = X.x+Y.x;
tmp.y = X.y+Y.y;
tmp.z = X.z+Y.z;
return(tmp);
}
vec3 scalar_vec3mult( double X, vec3 Y)
{
vec3 tmp;
tmp.x = X*Y.x;
tmp.y = X*Y.y;
tmp.z = X*Y.z;
return(tmp);
}