---
tags: synchronization
---
# Paul's M1 internship

###### tags: `synchronization`

---

## Background

A particle confined between two sinusoidally vibrated plates can either synchronize with the vibrations of the plate or stay in a chaotic state.

Here is a drawing of the plate:


The most important parameters are:

- $A$ the amplitude of the vibration (how much the box can move)
- $f$ the frequency of oscillation
- $h$ the height of the box
- $g$ the gravitational constants
- $\sigma$ the diameter of the particle
- The dissipative properties of the particle

______

For example, here, we see the z(t) of **non interacting** particles with time. They are confined between the oscillating plates drawn as the sinusoids:

| synchronized state | State with period doubling | chaotic state |
| --- | -------------------------------------------------- | -------------------------------------------------- |
| you can imagine it (i don't have any data ) |  |  |

The transition from the synchronized state to the asynchronized state depends on all the parameters listed above.

Here is for example a "synchronization map" (for additional information;, for example about the defiition of $s$,[ see this link](https://www.dropbox.com/scl/fi/37abkgvvfpdmae45s4rl0/absorbing-phase-transition-with-synchronization.pdf?rlkey=3f56y6tslphxemqvj6uu5hl4h&st=yq62bawg&dl=0)) of the system for various amplitude and height (at fixed $g$, $f$, ...):


We see that there is a chaotic attractor (asynchronized region $\equiv s < 1$) for intermediate amplitude.

If we do a cut at a fixed $h$ and look at the position of $N$ particles at the same given time, we obtain a chaotic figure close to a logistic map:


The stroboscopic position vs velocity phase space of one particle is also very nice!

## Plans

The goal for now, is to simulate and describe theoretically this phenomenon. The particle will be modeled as a hard disk which collides with oscillating plates following the usual collision rule:

$$v' = v + (1+ \alpha)(v_{plate} - v)$$

where $v'$ and $v$ are respectively the post and pre collisional velocity of the particle. $\alpha$ is a coefficient of restitution (proportional to the energy loss at collision due to dissipation). $v_{plate}$ is the velocity of the plates which follows from $v_{plate} = \dfrac{d x_{plate}(t)}{dt}$ and $x_{plate}^{bottom}=A\sin(2\pi f t)$, $x_{plate}^{top}=A\sin(2\pi f t) + h$

The particle should be affected by gravity too. The diameter of the particle can be get rid off by correct normalization of the height of the plate.

## To Do

::: info
- Code your own simulation of this phenomenon for $N$ **non-interacting** particles. Preferably in something fast (C, C++, rust, julia, python w/ numba...) so that you can easily iterate your results. Here is an (awful) time driven example code in rust (main directory) or in C (directory "C code"): https://github.com/Syrocco/synchroModel. You can also code a (faster) event/collision driven algorithm. Try to obtain a non chaotic -> chaotic -> non chaotic behavior as given in the figure above by varying $A$ for example.
:::

::: info
- This done, you should want to understand how and when synchronization arises as a function of the parameters of your system. To do this, you will want to analyze the collisional map of the system which gives you the time of the next collision knowing the time of the last one (and position of the collision). It is known that a bouncing ball (wihtout roof) exhibits chaotic behavior. [**See this link**](https://www.dropbox.com/scl/fo/51urpmin7z4yyqqixqt76/APUFWp_pKOWWtAkwjoFqiP8?rlkey=6drlp15jo2y06iqza140smg7r&st=9uvt6cjc&dl=0) for ressources on the theoretical analysis of the stability of the system see also these (**[1](https://vrs.amsi.org.au/wp-content/uploads/sites/92/2020/01/ceddia_julian_vrs-report.pdf), [2](https://sci-hub.ru/10.1103/physreve.48.3988 ), [3](https://iopscience.iop.org/article/10.1088/1742-6596/246/1/012003/pdf)**) references . You will want to do the same analysis done in the pdf of the dropbox, but for your system you might want to read [**this reference first**](https://www.dropbox.com/scl/fi/7olusggpaao46fbib391w/stability_maps.pdf?rlkey=7u6sxhxp7qj9h1jhli1fiumca&st=hkrxbqw7&dl=0), which is in 1D and very easy to follow. The complete system is probably too hard to analyze. You can start with a vanishing $g$. This will make the analysis easier (and it looks like that we can still observe the chaotic dynamics). You might also want to assume that $A$ is small (but $f$ large) so that the position of the plate does not change much and perhaps that $\alpha = 1$. In these conditions, the collision map is given by:
$$v_{0\rightarrow 1} = v_0~,~t_0 = 0$$
$$v_{1\rightarrow 2} = v_{0\rightarrow 1} - (1+\alpha)(v_{0\rightarrow 1} - A\omega\cos(\omega t_1))~,~t_1 = t_0 +|v_{0\rightarrow 1}|/z_0$$
$$v_{2\rightarrow 3} = v_{1\rightarrow 2} - (1+\alpha)(v_{1\rightarrow 2} - A\omega\cos(\omega t_2))~,~t_2 = t_1 + |v_{1\rightarrow 2}|/h$$
$$v_{3\rightarrow 4} = v_{2\rightarrow 3} - (1+\alpha)(v_{2\rightarrow 3} - A\omega\cos(\omega t_3))~,~t_3 = t_2 + |v_{2\rightarrow 3}|/h$$
$$...$$
$$v_{n\rightarrow n+1} = -\alpha v_{n-1\rightarrow n} + A\omega (1+\alpha)\cos(\omega t_n)~,~t_n = t_{n - 1} + |v_{n-1\rightarrow n}|/h$$

This is very reminiscent of the standard map.
:::

ideas:
#include<stdio.h>
#include <math.h>
#include <stdlib.h>
#include<time.h>
#define M_PI 3.1415926535
double vcol1(double A,double f,double x0,double v0,double t,double t0,double res);
double initp(double h);
double pseudorand(double max);
double trajx(double x0,double v0,double t);
double trajv(double x0,double v0,double t);
int binarySearch(double arr[], int left, int right, double target, double precision);
int main(){
double h=0.1;
int N=2000;
double dt=0.00001;
int i,k;
double A=h/100;
double f=2500;
double res =1;

double x[N];
x[0]=0.05;
double v[N];
v[0]=0;
double t[N];
t[0]=0;
printf("%lf " ,x[0]);
for (i=1;i<N;i++){
double Dr1=pow(v[i-1]/9.81,2)-(2/9.81)*(h-x[i-1]);
if ((Dr1>=0) && (v[i-1]>0)){
if (v[i-1]/9.81-pow(Dr1,0.5)>0){
t[i]=t[i-1]+v[i-1]/9.81-pow(Dr1,0.5);
x[i]=h-h/2000;}
else {
double Dr2=pow(v[i-1]/9.81,2)+(2/9.81)*(x[i-1]);
if (Dr2>0){
t[i]=t[i-1]+v[i-1]/9.81+pow(Dr2,0.5);
v[i]=vcol1(A,f,x[i-1],v[i-1],t[i],t[i-1],res);
// if (v[i]<0){v[i]=0.001;}
x[i]=h/2000;}

else {
printf("ERREUR -1");
exit(3);
return -1;
}}

}
else {
double Dr2=pow(v[i-1]/9.81,2)+(2/9.81)*(x[i-1]);
if (Dr2>=0){
t[i]=t[i-1]+v[i-1]/9.81+pow(Dr2,0.5);
v[i]=vcol1(A,f,x[i-1],v[i-1],t[i],t[i-1],res);
//if (v[i]<0){v[i]=0.001;}
x[i]=h/2000;
}
else {
printf("ERREUR -1");
exit(3);
return -1;
}

}

}
FILE* fichier = NULL;
double xg[10*N];
double tg[10*N];
fichier = fopen("val.txt", "w");
for (i=0;i<N;i++){
xg[10*i]=x[i];
tg[10*i]=t[i];
fprintf(fichier, "%lf %lf\n", xg[10*i],tg[10*i]);
dt=(t[i+1]-t[i])/9;
for (k=1;k<10;k++){
xg[10*i+k]=trajx(x[i],v[i],k*dt);
tg[10*i+k]=t[i]+k*dt;
fprintf(fichier, "%lf %lf\n", xg[10*i+k],tg[10*i+k]);
}

}
fclose(fichier);

return 1;
}
double vcol1(double A,double f,double x0,double v0,double t,double t0,double res){
double v=trajv(x0,v0,t-t0);
return v+(1+res)*(2*M_PI*f*A*cos(2*M_PI*f*t)-v);
}

double trajx(double x0,double v0,double t){
return (-9.81/2*pow(t,2)+x0+v0*t);}
double trajv(double x0,double v0,double t){
return v0-9.81*t;
}
double initp(double h){
return pseudorand(h);
}
double pseudorand(double max)
{
srand((unsigned) time(NULL));
return (max / RAND_MAX) * rand();}

first try to obtain a good simulation:
problems:
1:need to place the sphere at an arbitrary distance of the boundary otherwise it buggs "inside" the boudary and gives like 2 rebounds before crashing
2: this arbitrary distance does affect the simulation by a stupidly crazy factor, a change of ~1/1000 with otherwise the same parameters giving 2 very different simulation
d=h/1000

d=h/2000

next: test with a triangular signal for the plates

here is a bifurcation diagramm for A
It is not continuous but i think the triangular drive can cause this

we see with the following that, surprisingly, f will not play any role.(here a try with A=0.25 and h=1)
now a test for h with A=0.15 and h between 0.15 and 20
not very clear here but it seems like it alternate between bands of xhaos and periodicity
maybe a zoom
but i am not to excited to call it chaos as it seems to be more of an alternation between little number cycles to big ones (for instance arround 18.5 we go from a 5-cycle to a 12-ish-cycle) the big ones are nontheless far less stable that the ones we see in the non chaotic region so it is certainly as far as we can get with this.

here the restitution coefficient is 1
which gives me the idea of a diagramm with this as a parameter

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

double precision = 0.00000001;

int N = 10000;

double start(double h, double A, double f, int p, double x, double v);
double inter(double h, double A, double f, double x, double v, double start, int p);
double vcol(double res, double A, double f, double v0, double t);
double pos(double h, double A, double f, int p, double t);
void process_data(int Ng, double f, double x[], double tg[], int REB, double xt[], double Y[], double PARA[]);

int main() {
int Ng = 10000;
int REB = 50;
double *xg = malloc(Ng * 1000 * sizeof(double));
double *tg = malloc(Ng * 1000 * sizeof(double));
double *PARA = malloc(Ng * 1000 * sizeof(double));
double *xt = malloc(Ng * REB * sizeof(double));
double *Y = malloc(Ng * REB * sizeof(double));
double x[N];
double v[N];
double t[N];
double tb[N];

v[0] = 1; // must be non-zero
t[0] = 0;
tb[0] = 0;
v[0] = 2;
double h = 0.1;
double f = 5;
double res = 0.75;

int k = 0;
for (int l = 1; l <= Ng; l++) {
double A = 0.0001 * l;
double para = A;
x[0] = -A; // must be out of [A, h - A]
for (int i = 0; i < N - 1; i++) {
int p = -1;
if (x[i] == pos(h, A, f, 0, t[i])) {
if ((x[i] < A) && (x[i] < h - A)) {
if (v[i] < 0) {
p = 0;
}
if (v[i] > 0) {
if (((A - x[i]) / v[i] < (floor(t[i] * f) + 1.5) / f - t[i]) || ((h - A - x[i]) / v[i] < (floor(t[i] * f) + 1) / f - t[i])) {
p = 1;
} else {
p = 0;
}
}
}
if ((x[i] < A) && (x[i] > h - A)) {
if (v[i] > 0) {
p = 1;
}
if (v[i] < 0) {
if ((h - A - x[i]) / v[i] < (floor(t[i] * f) + 0.5) / f - t[i]) {
p = 0;
} else {
p = 1;
}
}
}
}
if (x[i] == pos(h, A, f, 1, t[i])) {
if ((x[i] > h - A) && (x[i] > A)) {
if (v[i] > 0) {
p = 1;
}
if (v[i] < 0) {
if (((h - A - x[i]) / v[i] < (floor(t[i] * f) + 1) / f - t[i]) || ((A - x[i]) / v[i] < (floor(t[i] * f) + 0.5) / f - t[i])) {
p = 0;
} else {
p = 1;
}
}
}
if ((x[i] > h - A) && (x[i] < A)) {
if (v[i] < 0) {
p = 0;
}
if (v[i] > 0) {
if ((A - x[i]) / v[i] < (floor(t[i] * f) + 0.5) / f - t[i]) {
p = 1;
} else {
p = 0;
}
}
}
}
double s = start(h, A, f, p, x[i], v[i]);
t[i + 1] = inter(h, A, f, x[i], v[i], s + t[i], p);
tb[i + 1] = (t[i + 1] * f - floor(t[i + 1] * f));
v[i + 1] = vcol(res, A, f, v[i], t[i + 1]);
x[i + 1] = pos(h, A, f, p, t[i + 1]);

if (i >= N - 1001) {
xg[k] = x[i];
tg[k] = t[i];
PARA[k] = para;
k++;
}
}
}

process_data(Ng, f, xg, tg, REB, xt, Y, PARA);

FILE *file = fopen("val.txt", "w");
if (file == NULL) {
perror("Unable to open file");
free(xg);
free(tg);
free(PARA);
free(xt);
free(Y);
return 1;
}

for (int i = 0; i < Ng * REB; i++) {
fprintf(file, "%.15lf %.15lf\n", xt[i], Y[i]);
}

fclose(file);
free(xg);
free(tg);
free(PARA);
free(xt);
free(Y);
return 0;
}

double start(double h, double A, double f, int p, double x, double v) {
if (x < A) {
if (p == 1) {
return (h - A - x) / v;
} else {
return 0;
}
}
if (x > h - A) {
if (p == 0) {
return (A - x) / v;
} else {
return 0;
}
}
return 0;
}

double inter(double h, double A, double f, double x, double v, double start, int p) {
double a = start;
int d = floor(a * f * 2);
double b = d + 1 / (2 * f);
double t = -1;

if (p == 0) {
if (d % 2 == 0) {
if (x < A) {
t = (x + A * (1 + 2 * d) - v * a) / (4 * A * f - v);
} else {
t = (A * (2 + 2 * d) - v * a) / (4 * A * f - v);
}
} else {
if (x < A) {
t = (x + A * (1 + 2 * (d + 1)) - v * a) / (4 * A * f - v);
} else {
if ((-2 * A) / v <= (d + 1) / (2 * f) - a) {
t = (-A * (2 * d) - v * a) / (-4 * A * f - v);
} else {
t = (A * (2 + 2 * d) - v * a) / (4 * A * f - v);
}
}
}
}
if (p == 1) {
if (d % 2 == 1) {
if (x > h - A) {
t = (x - h - A * (1 + 2 * d) - v * a) / (-4 * A * f - v);
} else {
t = (-A * (2 + 2 * d) - v * a) / (-4 * A * f - v);
}
} else {
if (x > h - A) {
t = (x - h - A * (1 + 2 * (d + 1)) - v * a) / (-4 * A * f - v);
} else {
if ((2 * A) / v <= (d + 1) / (2 * f) - a) {
t = (A * (2 * d) - v * a) / (4 * A * f - v);
} else {
t = (-A * (2 + 2 * (d + 1)) - v * a) / (-4 * A * f - v);
}
}
}
}

return t;
}

double vcol(double res, double A, double f, double v0, double t) {
int i = floor(t * f * 2);
if (i % 2 == 0) {
return v0 + (1 + res) * (4 * A * f - v0);
} else {
return v0 + (1 + res) * (-4 * A * f - v0);
}
}

double pos(double h, double A, double f, int p, double t) {
int i = floor(t * f * 2);
if (p == 0) {
if (i % 2 == 0) {
return 4 * A * f * (t - (i) / (f * 2)) - A;
} else {
return -4 * A * f * (t - (i) / (f * 2)) + A;
}
} else {
if (i % 2 == 0) {
return 4 * A * f * (t - (i) / (f * 2)) - A + h;
} else {
return -4 * A * f * (t - (i) / (f * 2)) + A + h;
}
}
}

void process_data(int Ng, double f, double x[], double tg[], int REB, double xt[], double Y[], double PARA[]) {
int j, i, k;
double b;

for (j = 1; j < Ng; j++) {
b = floor(tg[1000 * j - 1] * f);
k = 0;

if (b < 0) {
j++;
if (j >= Ng) break; // Ensure j doesn't go out of bounds
b = floor(tg[1000 * j - 1] * f);
}

for (i = 0; i < REB; i++) {
while (tg[1000 * j - 1 - k] > (b - i) / f) {
k++;
if (1000 * j - 1 - k < 0) break; // Ensure k doesn't go out of bounds
}

int idx1 = 1000 * j - 1 - k;
int idx2 = 1000 * j - 1 - (k - 1);

if (idx1 < 0 || idx2 < 0) break; // Ensure indices are valid

xt[REB * (j - 1) + i] = x[idx1] +
((x[idx2] - x[idx1]) /
(tg[idx2] - tg[idx1]) *
((b - i) / f - tg[idx1]));

Y[REB * (j - 1) + i] = PARA[idx1];
}
}
}

::: success
Raphaël's excursion:

Triangular signal, no gravity, time driven:


strong dependence of chaotic behavior on $\alpha$. And almost no dependence on full chaos or full synchro on $dt$, but strong dependence on edge of chaos close to the period doubling/quadrupling/etc region
:::
![Uploading file..._nzivwi78c]()



here a simulation with h=1 A=1 and res=0.5

here a simulation with h=0.5 A=0.5 and res=0.99
A>h seems to break the simulation down
like this with A=0.6
it comes down to the way i coded it , i'll try to improve that

problem solved

code qui marche pas

# -*- coding: utf-8 -*-
"""
Created on Mon Jun 3 10:17:23 2024

@author: 33695
"""

import pandas as pd

import xlrd
from numpy import loadtxt
import math
import numpy as np
import matplotlib.pyplot as plt
import random
import statistics
N=1000
f= loadtxt("traj.txt")
"""
print(f)
xg=[-1]*N*10
tg=[-1]*N*10

for i in range (5*N):
xg[i]=f[2*i]
tg[i]=f[2*i+1]
print(xg[i],tg[i])"""
F=1
A=1
h=1
N=10
Nt=9

t = f[0:N,1]
x=f[0:N,0]
v=f[0:N,2]
tg=[-1]*((len(t)-1)*10+1)
xg=[-1]*((len(t)-1)*10+1)
def para(x0,v0,a,A,F,h):
return v0 * (a) + x0
trb=[0]*2*Nt
trh=[0]*2*Nt
ttri=[0]*2*Nt

for i in range(Nt):
trb[2*i]=-A
trb[2*i+1]=A
trh[2*i]=h-A
trh[2*i+1]=h+A
ttri[2*i]=2*i/(2*F)
ttri[2*i+1]=(2*i+1)/(2*F)

plt.plot(ttri,trb)
#plt.plot(ttri,trh)
for i in range (len(t)):
tg[10*i]=t[i]
xg[10*i]=x[i]
if i<len(t)-1:
for k in range(1,10):
tg[10*i+k]=t[i]+k*(t[i+1]-t[i])/10
xg[10*i+k]=x[i]+v[i]*(tg[10*i+k]-t[i])-9.81/2*(tg[10*i+k]-t[i])**2

plt.plot(tg,xg)
#plt.plot(tp,xp)
"""trb=[0]*2*N
trh=[0]*2*N
ttri=[0]*2*N
for i in range(N):
trb[2*i]=-A
trb[2*i+1]=A
trh[2*i]=h-A
trh[2*i+1]=h+A
ttri[2*i]=2*i/(2*F)
ttri[2*i+1]=(2*i+1)/(2*F)
plt.plot(ttri,trb)
plt.plot(ttri,trh)"""
#plt.xlim(1495,1499)
#plt.scatter(fr,tg,s=3,alpha=0.5)
plt.show()





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

const double NEVER = 100000000000;
const double PRECISION = 0.00000001;
double g = 9.81;
int N = 10;

// Function prototypes
double start(double h, double A, double f, int p, double x, double v);
double lim(double x, double v, double t, double f, double A, double h, int p);
double vcol(double res, double A, double f, double v0, double t, double t0);
double pos(double h, double A, double f, int p, double t);
double fleche(double x0, double v0);
double seg(double A, double h, double f, double x, double t, double v, int p);
double collage(double A, double f, int p, double t);
void para(double t0, double v, double x, double A, double tp, double f, double p, double h, int i, double *r1, double *r2);

int main() {
double x[N];
double v[N];
double res = 0.6;
double t[N];
t[0] = 0;
x[0] = 1;
double h = 5461698651948616;
double A = 1;
double f = 1;

FILE *file = fopen("sin.txt", "w");
if (file == NULL) {
perror("Unable to open file");
return 1;
}
fprintf(file, "%.15lf %.15lf %.15lf \n", x[0], t[0], v[0]);

for (int i = 0; i < N - 1; i++) {
double th = NEVER;
int p = 0;
double tb = seg(A, h, f, x[i], t[i], v[i], 0);

if (fleche(x[i], v[i]) > h - A) {
th = seg(A, h, f, x[i], t[i], v[i], 1);
}

if ((tb < th) && (tb > t[i])) {
t[i + 1] = tb;
p = 0;
} else {
if ((th < tb) && (th > t[i])) {
t[i + 1] = th;
p = 1;
}
}

v[i + 1] = vcol(res, A, f, v[i], t[i + 1], t[i]);
x[i + 1] = pos(h, A, f, p, t[i + 1]);
/*
int test = floor(t[i] * 2 * f);
if ((test % 2 == 1) && (v[i] <= 16 * A * (2 * (t[i] - test / (2 * f)) - 1 / (2 * f)) * (1.0001)) && (x[i] < pos(h, A, f, 0, t[i]) + A / 100)) {
t[i + 1] = collage(A, f, 0, t[i]);
x[i + 1] = pos(h, A, f, 0, t[i + 1]);
v[i + 1] = 4 * A * f;
}
if ((test % 2 == 0) && (v[i] <= -16 * A * (2 * (t[i] - test / (2 * f)) - 1 / (2 * f)) * (1.0001)) && (x[i] < pos(h, A, f, 0, t[i]) + A / 100)) {
t[i + 1] = collage(A, f, 0, t[i]);
x[i + 1] = pos(h, A, f, 0, t[i + 1]);
v[i + 1] = 4 * A * f;
}
if (t[i + 1] < t[i]) {
t[i + 1] = collage(A, f, 0, t[i]);
x[i + 1] = pos(h, A, f, p, t[i + 1]);
}*/

// Output the results to file and console
fprintf(file, "%.15lf %.15lf %.15lf \n", x[i + 1], t[i + 1],v[i + 1]);
printf("xg %lf v %lf tg %lf \n", x[i + 1], v[i + 1], t[i + 1]);
}

fclose(file);
return 0;
}

// Velocity after collision
double vcol(double res, double A, double f, double v0, double tnow, double tlastcoll) {
int i = floor(tnow * f * 2);
double v = v0 - g * (tnow - tlastcoll);
double t0 = i / (2 * f);
int sign = -1;
if (i % 2 == 0) {
sign = 1;
}
double vplate = sign*8*A*f*(1 + 4*f*(t0 - tnow));
return v + (1 + res) * (vplate - v);
}

// Position function
double pos(double h, double A, double f, int p, double t) {
int i = floor(t * f * 2);
if (p == 0) {
if (i % 2 == 0) {
t=t-i/(2*f);
return -16 * f * f * A * (t ) * (t - 1 / (2 * f));
} else {
t=t-(i+1)/(2*f);
return 16 * f * f * A * (t - i / (2 * f)) * (t - i / (2 * f) - 1 / (2 * f));
}
} else {
if (i % 2 == 0) {
t=t-i/(2*f);
return -16 * f * f * A * (t - i / (2 * f)) * (t - i / (2 * f) - 1 / (2 * f)) + h;
} else {
t=t-(i+1)/(2*f);
return 16 * f * f * A * (t - i / (2 * f)) * (t - i / (2 * f) - 1 / (2 * f)) + h;
}
}
}

// Compute the height reached by the ball
double fleche(double x0, double v0) {
return v0 * v0 / (g * 2) + x0;
}

double seg(double A, double h, double f, double x, double t, double v, int p) {
int n = 0;
//while (n / (2 * f) < 500) {
while (1){
double t0 = (floor(t * 2 * f) + n) / (2 * f);
int a = (floor(t * f * 2) + n);

int sign = 1;
if ((a % 2) == 0){
sign = -1;
}

double o = -g / 2.0 - 16.0 * A * f * f * sign;
double b = g * t + v + 32.0 * A * f * f * t0 * sign + 8.0 * A * f * sign;
double c = -g / 2.0 * t * t - v * t + x - 16.0 * A * f * f * t0 * t0 * sign - 8.0 * A * f * t0 * sign + p * h;
double d = b * b - 4 * o * c;
printf("n: %d, t0: %.2lf, D: %lf\n", n, t0, d);
if (d >= 0) {
double t1 = (-b + sqrt(d)) / (2 * o);
double t2 = (-b - sqrt(d)) / (2 * o);

if ((t2 > t0) && (t2 <= t0 + 1 / (2 * f)) && (t2 > t + 1 / (100/PRECISION * f))) {
return t2;
} else {
if ((t1 > t0) && (t1 <= t0 + 1 / (2 * f)) && (t1 > t + 1 / (100/PRECISION * f))) {
return t1;
} else {
n += 1;
}
}
} else {
n += 1;
}

}

}

// Calculate the time limit
double lim(double x, double v, double t, double f, double A, double h, int p) {
return ((v + g * t) / g + sqrt((v + g * t) * (v + g * t) + 2 * g * (x - v * t + A - p * h - g * t * t / 2)) / g);
}

double collage(double A, double f, int p, double t) {
int i = floor(t * f * 2);
if (i % 2 == 0) {
return (i + 1 + PRECISION) / (2 * f) - 1 / (4 * f);
} else {
return (i + 2 + PRECISION) / (2 * f) - 1 / (4 * f);
}
}


differentes frequences
jolie hyperbole
almost static plates.

moi je trouve que la divergence est pas mal quand meme

belle
le plateau au debut c'est la fin de la simulation