#include <stdio.h>
#include <math.h>
#include "rtbp.h"
#include "ctants.h"
#include "nwtnev.h"

struct proptraj_fdf_prm {
    double r0, v0;
    int nt;
    double pmax;
    double dm;
};

void proptraj_fdf (double alf, double *u, double *dalfu, void *prm) {
    struct proptraj_fdf_prm *p = prm;
    double r0=(*p).r0, v0=(*p).v0, nt=(*p).nt, pmax=(*p).pmax;
    
    double z[4], zf[4], dz[4];
    rva2z(MU, r0, v0, cos(alf), sin(alf), z);
    proptraj(MU, nt, z, zf, dz, pmax, NULL);
    double rzf = (MU - 1.0) - zf[0];
    double dm = rzf * AMOON;
    if (dm < 0.0) dm += RMOON;
    else dm -= RMOON;
    (*p).dm = dm;
    *u = zf[2];
    *dalfu = dz[0]*(-r0*sin(alf)) + dz[1]*(r0*cos(alf)) + dz[2]*(-v0*cos(alf)) + dz[3]*(-v0*sin(alf));
}

int main (int argc, char *argv[]) {
    double tolu, pmax;
    int maxit, ixrr;
    if (argc<5
         || sscanf(argv[1], "%lf", &tolu)!=1
         || sscanf(argv[2], "%d", &maxit)!=1
         || sscanf(argv[3], "%lf", &pmax)!=1
         || sscanf(argv[4], "%d", &ixrr)!=1
      ) {
        fprintf(stderr,"%s tolu maxit pmax ixrr\n", argv[0]);
        return -1;
    }

    double r0, v0, alf0, nt;
    while (fscanf(stdin, "%lf %lf %lf %lf", &r0, &v0, &alf0, &nt)==4) {
        struct proptraj_fdf_prm prm;
        prm.r0=r0;
        prm.v0=v0;
        prm.nt=nt;
        prm.pmax=pmax;
        if (newton(&alf0, tolu, -1, maxit, proptraj_fdf, &prm, ixrr) != -1) {
            fprintf(stdout, "%.16g %.16g %.16g %.16g\n", r0, v0, alf0, prm.dm);
        } else if (ixrr) {
            fprintf(stderr, "No s'ha trobat cap solució\n");
        }
    }
    
    return 0;
}