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

#define pi              M_PI
//#define nombre1         "1D1c.dat"

/*  El valor de LARGO debe ajustarse para asegurar convergencia.     */
/*  Con eta cerca de 1.0  LARGO debe ser bastante grande             */

#define    LARGO          3000
#define    grav           1.0            /* acel de gravedad            */
#define    AA             1.0            /* amplit del suelo            */
#define    ene            200

FILE       *archivo1;

short int  abortar,hayraiz;
int        nn,kk,xj;
long int   cont, contS,contR;
double     Ohm,yy1,v1,z,vz,ts,tsim,tsemilla,tfinal,Ts,nu,eta;
char	   nombre[30];
/***********************************************************************/
/*  Esta funciOn se anula cuando la partIcula toca el suelo           */
double fz(double tx)
{ double taux,zz;
    taux = tsim + tx;
    zz   = AA*sin(Ohm*taux);
    return(yy1 + v1*tx - 0.5*grav*tx*tx - zz);
}
void Newton()
{ double t0, t1, f0, f1, Af1;
  int mucho;
    hayraiz = 0;
    t0      = tsemilla;
    f0      = fz(t0);
    t1      = t0 + 1e-3/Ohm;
    f1      = fz(t1);
    mucho   = 0;
    do
    {if(f1!=f0) tfinal = (t0*f1-t1*f0)/(f1-f0);
     t0     = t1;
     f0     = f1;
     t1     = tfinal;
     f1     = fz(t1);
     Af1    = fabs(f1);
     mucho++;
    }while(  (Af1>1e-8) && (mucho<50)  ) ;
    if(mucho<50) hayraiz=1;
}
/***********************************************************************/
void inicial()
{int jj;
    abortar  = 0;
    tsim     = pi/Ohm;
    vz       = AA*Ohm*cos(Ohm*tsim);
    yy1      = AA*sin(Ohm*tsim) + 0.1e-6;
    v1       = vz + 20.0*drand48()*grav/Ohm;
}
/***********************************************************************/
void calculo_ts()
{ double taux,tau,yaux,zaux;
    abortar= 0;
    tau    = 0.1/Ohm;      /*     "peque~no" salto de tiempo             */
    taux   = 0.0;
    do
    {do                    /* loop que avanza hasta quedar bajo el suelo */
     {taux   = taux + tau;
      zaux   = AA*sin(Ohm*(tsim+taux));
      yaux   = yy1 + v1*taux - 0.5*grav*taux*taux;
     }while(yaux>zaux);
     tsemilla  = taux;
     Newton();
     ts        = tfinal;

     tau = 0.7*tau;
     if(ts<1e-5)  contS++;
     else contS=0;             /* contS = veces que bote es muy breve*/
     if(  (contS>20) || (hayraiz==0) )
     { abortar=1;
       printf("abortado:    contS=%4d  hayraiz=%d \n",contS,hayraiz);
       break;
     }
    }while(ts<0 && abortar==0);
}
/***********************************************************************/
void suelo()
{tsim  = tsim + ts;
 z     = AA*sin(Ohm*tsim);
 vz    = AA*Ohm*cos(Ohm*tsim);
 yy1   = z;
 v1    = v1 - ts*grav;
 v1    = (1+eta)*vz - eta*v1;
 cont++;
}
/***********************************************************************/
void Evolucion()
{ cont      = 0;
  contS     = 0;
  inicial();
  while(cont<=LARGO)
  { calculo_ts();
    if(abortar) break;
    suelo();
  }
}
/************************************************************************/
void postEvolucion()
{ double  v10,t0,ttot,nu;
  int     jj;
  v10  =  v1;
  t0   =  tsim;
  while(cont<=LARGO+ene)
  { calculo_ts();
    if(abortar) break;
    suelo();
    nu = (tsim-t0)/Ts;
    if(nu>1.8)
    { archivo1=fopen(nombre,"at");
      contR++;
         fprintf(archivo1,"%7.4f %8.4f\n",Ohm,v1);
      fclose(archivo1);
    }
  }
}
/***********************************************************************/
main()
{nn    = 0;
 contR = 0;
 srand48(100);
 eta = 0.15; 
 do
 {  Ohm = 0.9;
    sprintf(nombre,"birfur%g",eta); //se quedar'a pa siempre con "bir"
    archivo1=fopen(nombre,"wt");
    fprintf(archivo1," # Ohm  v1\n");
    fclose(archivo1);
    do
    {  abortar    = 0;
       Ohm       += 0.00025;
       Ts         = 2*pi/Ohm;
       Evolucion();
       postEvolucion();
       if(nn%50==0) printf("W=%f\n",Ohm);
       nn++;
    }while(Ohm<1.8);
    eta += 0.05;
 }while(eta<0.96);
}