// C++ Program

#include <iostream>
#include <math.h>       /* pow */
#include <fstream>

using namespace std;
// For real declaration 

typedef double R;
const double pi = 3.141592654 ; 
const double g = 9.81 ;
const double rho = 2700. ;
const double v0 = 0.0 ;
const int nb_it = 200000000 ;
const int freq_sauv = 5000 ;
//
const int Np = 1 ;
const double h0[Np] = {1.} ;
const double Ri[Np] = {0.15} ;
const double x0[Np] = {0.75} ;
const double rest[Np] = {0.7} ;
//
const double Young_p = 7.0e+10 ; // Young's modulus of the particle 
const double Young_b = 7.0e+10 ; // Young's modulus of the wall 
const double Poiss_p = 0.3 ; // Poisson's ratio of the rigid wall
const double Poiss_b = 0.3 ; // Poisson's ratio of the spherical particle
const double Inv_E = (1.- pow(Poiss_p,2))/Young_p + (1.- pow(Poiss_b,2))/Young_b ; // Equivalent elastic coefficient particle/wall
const double E_eq = 1./Inv_E ; // Equivalent elastic coefficient particle/wall
//
const double timp = 373.15 ;
const double lambda_w = 237 ; // W/m.K Thermal conductivity of the material of the wall (aluminium) 
const double lambda_f = 0.025 ; // W/m.K Thermal conductivity of gas at 293 K (20°C) (air)
//
const double t0 = 298.15 ;
const double cap = 897 ; // J/kg.K Thermal capacity of the material of the ball (aluminium)
const double lambda = 237 ; // W/m.K Thermal conductivity of the material of the ball (aluminium)
//const double cap = 129 ; // J/kg.K Thermal capacity of the material of the ball (or)
//const double lambda = 318 ; // W/m.K Thermal conductivity of the material of the ball (or)
//const double cap = 377 ; // J/kg.K Thermal capacity of the material of the ball (bronze)
//const double lambda = 188.3 ; // W/m.K Thermal conductivity of the material of the ball (bronze)
//const double cap = 900 ; // J/kg.K Thermal capacity of the material of the ball (Céramique)
//const double lambda = 35 ; // W/m.K Thermal conductivity of the material of the ball (Céramique)
//
const double K_p = (1.- pow(Poiss_p,2))/Young_p ; // Equivalent elastic coefficient of the wall
const double K_b = (1.- pow(Poiss_b,2))/Young_b ; // Equivalent elastic coefficient of the particle
const double R_len = 1.09 ; // Thickness fluid ratio (fluid layer thickness surrounding the sphere) 
const double S_bar = 0.002 ; // Thickness fluid ratio (fluid layer thicknes between sphere and wall)
const double Delta_bar = 0e-2 ; // Thickness fluid ratio between particle and wall without contact



// FUNCTION Solution discrétisée
void Dep_Vit_Acc(double *zi , double *vi , double *zf , double *vf , double *ai , double *af , 
double *un, double *r, int it, double dt) {

double vn, m_b, psi, kn, fn ;

if (it == 0) {
	for (int i = 0 ; i < Np ; i++) {
	zi[i] = h0[i] ;
	vi[i] = v0 ;
	ai[i] = -g ;
	zf[i] = h0[i] ;
	vf[i] = v0 ;
	af[i] = -g ;
}
}
else {
	for (int i = 0 ; i < Np ; i++) {
	un[i] = zi[i]-Ri[i] ;
	if (un[i] < 0.) {
	vn = vi[i] ;
	m_b = rho*4*pi*pow(Ri[i],3)/3 ;
	psi = abs(log(rest[i]))/pow((log(rest[i])*log(rest[i])+pow(pi,2)),0.5) ;
	kn = E_eq*pow(Ri[i],0.5) ;
	fn = -4.*kn*pow(abs(un[i]),0.5)*un[i]/3. - 2.*pow(5./6.,0.5)*psi*pow(2.*kn*pow(abs(un[i]),0.5)*m_b,0.5)*vn ;
	af[i] = -g + fn/m_b ;
	r[i] = 3*(K_p+K_b)*Ri[i]*abs(fn)/4 ;
	r[i] = pow(r[i],0.33333333) ;
    }
    else {
    af[i] = -g ;
    r[i] = 0. ;
	}
	vf[i] = (af[i]+ai[i])*dt/2. + vi[i] ;
	zf[i] = af[i]*pow(dt,2)/2. + vf[i]*dt + zi[i] ;
	ai[i] = af[i] ;
	vi[i] = vf[i] ;
	zi[i] = zf[i] ;
}
}
}

// FUNCTION temperature à l'instant t
void Temperature_sphere(double *tf, double *ti , double &dt_th , double *r_cont , double *depth , double *fluhz , double *fluair , double *ratiohz , double *ratioair , int it, double dt) {
double conduc,delta,conduc_f,r_in,r_out,r_s, m_b ;

for (int i = 0 ; i < Np ; i++) {

if ( r_cont[i] == 0. || it == 0 ) {
	
conduc = 0. ;
conduc_f = 0. ;

	
	dt_th = dt ;
	tf[i] = t0 ;
	ti[i] = tf[i] ;
    fluhz[i] = 0. ;
    fluair[i] = 0. ;
    ratiohz[i] = 0. ;
    ratioair[i] = 0. ; 
}

else { 
	
	dt_th = dt*1000 ; 
	
	conduc = 4.*r_cont[i]/(1./lambda + 1./lambda_w) ;
	m_b = rho*4*pi*pow(Ri[i],3)/3 ;

	delta = -depth[i]/Ri[i] ;
    r_in = sqrt(1-pow((1-delta),2)) ;
    r_out = sqrt(pow(R_len,2)-pow((1-delta),2)) ;
    r_s = sqrt(1-pow((1-S_bar-delta),2)) ;
    conduc_f = (pow(r_s,2) - pow(r_in,2))/(2*S_bar) ;
    conduc_f = conduc_f + sqrt(1-pow(r_out,2)) ;
    conduc_f = conduc_f - sqrt(1-pow(r_s,2)) ;
    conduc_f = conduc_f + (1-delta)*log((1-delta-sqrt(1-pow(r_out,2)))/(1-delta-sqrt(1-pow(r_s,2)))) ;
     
    fluhz[i] = conduc*(timp-ti[i]) ;
    fluair[i] = 2.*pi*conduc_f*lambda_f*Ri[i]*(timp-ti[i]) ;

    ratiohz[i] = fluhz[i]/(fluhz[i] + fluair[i]) ;
    ratioair[i] = fluair[i]/(fluhz[i] + fluair[i]) ;
	
    tf[i] = ti[i] + dt_th*conduc*(timp-ti[i])/(cap*m_b) ;
    tf[i] = tf[i] + dt_th*2.*pi*conduc_f*lambda_f*Ri[i]*(timp-ti[i])/(cap*m_b) ;
    ti[i] = tf[i] ;
  } 

}
}

// FUNCTION critical time step
void Min_timestep(double &delta_min) {
double m_b[Np], psi, cn, cc, kn, delta_t ;
	
delta_min = 10 ;

for (int i = 0 ; i < Np ; i++) {
m_b[i] = rho*4*pi*pow(Ri[i],3)/3 ;
psi = abs(log(rest[i]))/pow((log(rest[i])*log(rest[i])+pow(pi,2)),0.5) ;
psi = abs(log(rest[i]))/pow((log(rest[i])*log(rest[i])+pow(pi,2)),0.5) ;
kn = E_eq*pow(Ri[i],0.5) ;
cn = pow((10*kn*m_b[i]/3),0.5)*psi ;
cc = 2*pow(kn*m_b[i],0.5) ;

delta_t = 2*sqrt(m_b[i]/kn)*(sqrt(1+pow(cn/cc,2))-cn/cc)/20. ;

if ( delta_t < delta_min ) {
	delta_min = delta_t ;
}
}

cout << "*****************************************" << endl ;
cout << "Pas de temps critique (calcul dynamique): " << delta_min << " s " <<endl ;
cout << "*****************************************" << endl ;

}

// MAIN PROGRAM
int main()
{
double z0[Np],v0[Np],zt[Np],vt[Np],a0[Np],af[Np],delta_t,tps,r_cont[Np] ;
double t_fin[Np],t_init[Np],depth[Np],flhz[Np],flair[Np],ratiohz[Np],ratioair[Np],dt_th ;
int    itt ;

Min_timestep(delta_t) ;

ofstream My_data0("chuteNp") ;
ofstream My_data1("chuteNp_at") ;
ofstream My_data2("tempNp") ;
ofstream My_data3("tempNp_at") ;
ofstream My_data4("champ_temp") ;
ofstream My_data5("flux_sph") ;

tps = 0. ;

for (itt = 0 ; itt <= nb_it ; ++itt) {
Dep_Vit_Acc(z0,v0,zt,vt,a0,af,depth,r_cont,itt,delta_t) ;
Temperature_sphere(t_fin,t_init,dt_th,r_cont,depth,flhz,flair,ratiohz,ratioair,itt,delta_t) ;

tps = tps + dt_th ;

if ( (itt%freq_sauv) == 0 ) {
My_data0 << itt/freq_sauv << " " ;
My_data2 << itt/freq_sauv << " " ;

for (int i = 0 ; i < Np ; i++) {	
My_data0 << x0[i] << " " ;
My_data2 << x0[i] << " " ;
}
for (int i = 0 ; i < Np ; i++) {	
My_data0 << z0[i] << " " ;
My_data2 << z0[i] << " " ;
}
for (int i = 0 ; i < Np ; i++) {	
My_data2 << t_fin[i]-273.15 << " " ;
}
for (int i = 0 ; i < Np ; i++) {	
My_data0 << Ri[i] << " " ;
My_data2 << Ri[i] << " " ;
}

My_data0 << endl ;
My_data0 << "  " << endl  ; 		
My_data0 << "  " << endl  ; 		
My_data2 << endl ;
My_data2 << "  " << endl  ; 		
My_data2 << "  " << endl  ; 	

My_data4 << tps/3600. << " " ;
for (int i = 0 ; i < Np ; i++) {	
My_data4 << t_fin[i]-273.15 << " " ;
}
for (int i = 0 ; i < Np ; i++) {	
My_data4 << flhz[i] + flair[i] << " " ;
}
for (int i = 0 ; i < Np ; i++) {	
My_data4 << flhz[i] << " " ;
}
for (int i = 0 ; i < Np ; i++) {	
My_data4 << flair[i] << " " ;
}
My_data4 << endl ;
	
}

if ( tps/3600 > 0.17 && tps/3600 < 0.18 ) {
for (int i = 0 ; i < Np ; i++) {	
My_data5 <<  i << " " << Ri[i] << " " << rho*4*pi*pow(Ri[i],3)/3. << " " << flhz[i] + flair[i] << " " << ratiohz[i] << " " << ratioair[i] << endl ;
}

}

}

My_data1 << tps << " " ;
My_data3 << tps << " " ;
for (int i = 0 ; i < Np ; i++) {	
My_data1 << x0[i] << " " ;
My_data3 << x0[i] << " " ;
}
for (int i = 0 ; i < Np ; i++) {	
My_data1 << z0[i] << " " ;
My_data3 << z0[i] << " " ;
}
for (int i = 0 ; i < Np ; i++) {	
My_data3 << t_fin[i]-273.15 << " " ;
}
for (int i = 0 ; i < Np ; i++) {	
My_data1 << Ri[i] << " " ;
My_data3 << Ri[i] << " " ;
}
My_data1 << endl ;
My_data3 << endl ;



}