//Solutions for Quiz 4, 2018

#include <cstdlib>
#include <cmath>
#include <cstdio>
#include <complex>
#include <vector>

const double PI=4.0*atan(1.0);
const double HBARC=197.3269602;

using namespace std;
double GetGamma(double Lratio){
	// Lratio=L*alpha/V0
	double Lcalc,dg,dLcalcdguess,root;
	double guess=pow(6.0*Lratio,1.0/3.0);
	do{
		root=sqrt(1.0+guess*guess);
		Lcalc=guess*root+asinh(guess);
		printf("gamma guess=%g, Lcalc-Lratio=%g\n",guess,Lcalc-Lratio);
		root=sqrt(1.0+guess*guess);
		dLcalcdguess=root+(guess*guess/root)+(1.0/root);
		guess=guess+(Lratio-Lcalc)/dLcalcdguess;
	}while(fabs(Lcalc-Lratio)>1.0E-10);
	Lcalc=guess*root+asinh(guess);
	printf("gamma=%17.15f, conv. test=%g =? 0\n",guess,Lcalc-Lratio);
	return guess;
}

void GetXYVXVY(double alpha,double L,double V0, double gamma,double t,double &x,double &y,double &vx,double &vy){
	double z,root;
	z=(gamma-alpha*t);
	root=sqrt(1.0+z*z);
	y=sqrt(1.0+gamma*gamma)*(V0/alpha)-root*(V0/alpha);
	vy=z*V0/root;
	vx=V0*sqrt(1.0+gamma*gamma)-V0*z*z/root;
	x=V0*sqrt(1.0+gamma*gamma)*t
		-(0.5*V0/alpha)*(asinh(z)-z*root)
			+(0.5*V0/alpha)*(asinh(gamma)-gamma*sqrt(1.0+gamma*gamma));
}




int main(int argc,char *argv[]){
	double V0=10.0,alpha=0.02,L=1000.0;
	double gamma,x,y,t,vx,vy;
	FILE *fptr=fopen("boat.dat","w");
	//gamma=0.8926677710351815; // works for L=100, v0=10, alpha=0.02
	gamma=GetGamma(L*alpha/V0);
	
	

	//for(gamma=0.892667771;gamma<0.89267773;gamma+=0.0000000001){
	
	// Calculate for tf
	t=2.0*gamma/alpha;
	printf("Max t=%g\n",t);
	GetXYVXVY(alpha,L,V0,gamma,t,x,y,vx,vy);
	printf("t=%12.7f, x,y=(%16.12f,%16.12f), vx,vy=(%g,%g)\n",t,x,y,vx,vy);
	//
	// Calculate for tf/2
	t=gamma/alpha;
	GetXYVXVY(alpha,L,V0,gamma,t,x,y,vx,vy);
	printf("t=%12.7f, x,y=(%16.12f,%16.12f), vx,vy=(%g,%g)\n",t,x,y,vx,vy);
	printf("Max y=%g\n",y);
	
	return 0;
}