// .............................................................. oscillforz.cc
// ................................................................... #include
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <grdriver.h>
#include "corso.h"

void func(double *t,double *y,double *dydt);

GTO to;
VideoData vd;

main()
{
  int DeltaSave,i,iy,iy0,iyf,iyf0,nPunti,nStep;
  double dt,tMax,tMin,yScale,yScale2;
  double *t,*tt,*v,*yStart;
  double **y;
  char answ[122],buff[40],InFile[80];
  FILE *hnd;

  tt=new double [8];
  t=new double [1000];
  v=new double [1000];
  yStart=new double [4];
  y=new (double *) [2];
  for (i=0;i<2;i++) y[i]=new double [1000];
  // .............................................................. costanti
  tt[1]=0.3; //        massa
  tt[2]=50.0; //       costante di Hooke
  tt[3]=0.4; //        smorzamento
  tt[4]=40.0; //       frequenza di guida
  tt[5]=5.0; //        ampiezza di guida
  // ................................................... condizioni iniziali
  yStart[0]=0.0; //    posizione iniziale
  yStart[1]=0.0; //    velocita' iniziale
  // ................................................... intervallo di tempo
  tMin=0.0;
  tMax=10.0;
  nStep=3200;
  // ......................................................... campionamento
  DeltaSave=4;
  dt=(tMax-tMin)*(double)(DeltaSave)/(double)(nStep);
  // ........................................................... Runge-Kutta
  RungKutCS(tt,yStart,2,tMin,tMax,nStep,t,y,DeltaSave,&nPunti,func);
  printf("nPunti=%d\n",nPunti);
  // ............................................................... grafico
  yScale=2000.0;
  yScale2=20.0;
  GrInstall();
  GrClearScreen(vd.brightwhite);
  GrHLine(0,799,300,vd.black);
  iy0=300-(int)(y[0][0]*yScale);
  iyf0=300;
  for (i=1;i<nPunti;i++)
    { iy=300-(int)(0.5+y[0][i]*yScale);
      iyf=300-(int)(0.5+tt[5]*sin(tt[4]*i*dt)*yScale2);
      GrLine(i-1,iy0,i,iy,vd.midred);
      GrLine(i-1,iyf0,i,iyf,vd.midblue);
      iy0=iy;
      iyf0=iyf;
    }
  fgets(answ,2,stdin);
  // ..................................................... libera la memoria
  for (i=1;i>=0;i--) delete [] y[i];
  delete [] y;
  delete [] yStart;
  delete [] v;
  delete [] t;
  delete [] tt;
}

// ---------------------------------------------------------------------------

void func(double *t,double *y,double *dydt)
{ 
  double a,k,m,omega,relax;

  m=t[1];
  k=t[2];
  relax=t[3];
  omega=t[4];
  a=t[5];

  // ................................................................. velocita'
  dydt[0]=y[1];
  // ............................................................. accelerazione
  dydt[1]=-(k/m)*y[0]-(relax/m)*y[1]+a*sin(omega*t[0]);
}