#define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include <math.h>

#define DT 20
#define T DT*1000
#define del 0.2
#define del2 0.2
#define omega1 sqrt(1-del*del)
#define omega2 sqrt(1-del2*del2)
#define Th 2.5
void main(void)
{
FILE *fp;
int i;
\define the initial conditions
double x = 0.1, t = 0.0, h = 1.0 / DT, y = 0.0,x0,y0;

fp = fopen("Kai.csv", "w+");
fprintf(fp, "tau, x, x'n");


for (i = 0; i <= T; i++) {
x0 = x; \ Save the new initial conditions for using when x' reach the boundary
y0 = y;
if (y > Th) {
    while (y > Th) {
        \Solution for x, x' with initial conditions x0, y0
        x = exp(-del*t)*((x0 - 1)*cos(omega1*t) + (1 / omega1)*(y0 + del*(x0 - 1))*sin(omega1*t)) + 1;
        y = exp(-del*t)*(y0*cos(omega1*t) - (1 / omega1)*(del*y0 + (x0 - 1))*sin(omega1*t));
        fprintf(fp, "%f, %f,%fn", t, x, y);
        t = t + h;
    }
}
else if (y < -Th) {
    while (y < -Th) {
        x = exp(-del*t)*((x0 + 1)*cos(omega1*t) + (1 / omega1)*(y0 + del*(x0 + 1))*sin(omega1*t)) - 1;
        y = exp(-del*t)*(y0*cos(omega1*t) - (1 / omega1)*(del*y0 + (x0 + 1))*sin(omega1*t));
        fprintf(fp, "%f, %f,%fn", t, x, y);
        t = t + h;
    }
}
else {
    while ((y>=-Th) && (y<=Th)) {
        x = exp(del2*t)*(x0*cos(omega2*t) + (1 / omega2)*(y0 - del2*x0)*sin(omega2*t));
        y = exp(del2*t)*(y0*cos(omega2*t) + (1 / omega2)*(del2*y0 - x0)*sin(omega2*t));
        fprintf(fp, "%f, %f,%fn", t, x, y);
        t = t + h;
        }
    }
}
fclose(fp);
}