Ball Throw Simulation

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

static const double PI      = 3.141592654;
static const double radDeg  = 57.2957795; // 180/PI Conversion factor rad to degrees

typedef struct vectorPolar
{
    double r;
    double theta;

} *vectorPolar;

void add(vectorPolar v1, vectorPolar v2)
{
    // Addition of two polar vectors r, theta.
    // Result stored in first vector.
    double x =  v1->r * cos(v1->theta) + v2->r * cos(v2->theta);
    double y =  v1->r * sin(v1->theta) + v2->r * sin(v2->theta);
    v1->r = sqrt(x * x + y * y);
    v1->theta = atan2(y, x);
}

void printDataPoint(int timeMin, int timeSeg, double height, vectorPolar p, vectorPolar v)
{
    // Copies all results of simulation to file 'coordenadas.txt'
    static FILE *filePointer = NULL;
    if (filePointer == NULL)
    {
        filePointer = fopen("coordenadas.txt", "w");
        fprintf(filePointer, "T(min), T(seg) , Height, Pos(magn), Pos(ang), Veloc(magn), Veloc(ang)\n");
    }
    fprintf(filePointer, "%d, %d, %E, %E, %E, %E, %E\n",
        timeMin, timeSeg, height, p->r, p->theta * radDeg, v->r, v->theta * radDeg);
}

int main (void)
{   // Simulation of a ball thrown from Earth surface
    // on a point at equator with given initial speed
    // and taking into account Earth's rotation
    // Simón Graffe simonjgl@gmail.com
    // constants:
    const double G              = 6.674E-11; // m^3.kg^-1.s^-2 Gravitational constant
    const double earthMass      = 5.9723E24; // kg             Earth's mass
    const double earthRadius    = 6378100;   // m              Equatorial radius
    const double earthRotPeriod = 0.99726968 * 3600 * 24; // d -> s  Rotation period
    const double GM = G * earthMass;         // m^3.s^-2       Precalculated factor
    const int    step           = 100;        // ms            Time step for simulation loop
    const double deltaT = (double)step/1000.; // s             Integration step

    struct vectorPolar p, v, a;

    // Initial position (on equator):
    p.r = earthRadius;
    p.theta = 0;

    // Tangent velocity at that point:
    struct vectorPolar vt;
    vt.r = 2 * PI * earthRadius / earthRotPeriod;
    vt.theta = PI / 2;

    // Throw speed (vertically)
    struct vectorPolar vi;
    vi.r = 38660E3/3600; // km/h -> m/s
    vi.theta = 0;

    // Combine both velocities:
    v.r = 0;
    v.theta = 0;
    add(&v, &vt); // v = v + vt
    add(&v, &vi); // v = v + vi -> v = vt +vi

    // Start of simulation:
    // Max time: 24 h in milli seconds
    for (int t = 0; t<=24*3600*1000 ;t += step)
    {
        // gravitational acceleration at that point
        a.r     = - GM / (p.r * p.r);
        a.theta = p.theta;
        struct vectorPolar deltaV;
        deltaV.r     = a.r * deltaT;
        deltaV.theta = a.theta;
        struct vectorPolar deltaP;
        deltaP.r     = v.r * deltaT;
        deltaP.theta = v.theta;
        // Update position and velocity at the new point
        add(&p, &deltaP);
        add(&v, &deltaV);

        double height; // Height above Earth surface
        height = (p.r - earthRadius);
        if (((t % 300000) == 0) || (height <0))
            printDataPoint((t/1000)/60, (t/1000) % 60, height, &p,&v);
        // If height is negative, the object crashed, we exit the simulation
        if (height < 0)
            break;
    }
    return 0;
}