#include<stdio.h>
#include<stdlib.h>
#include<math.h>
int main(void){
double omega;
double k;
int i;
double j;
double theta; /* Our theta is theta_inverse */
double alpha;
double denom1,denom2;
double f;
double tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8;
double c1, c2, c3, c4, c5, c6;
double Gp,nu,h,gamma,epsilon;
double delta, Em, Ep;

delta = 2.;
printf("# delta %le\n",delta);
gamma = 1./3.;
printf("# gamma %le\n",gamma);
nu = 1./3.;
printf("# nu %le\n",nu);
h = 11.061;
printf("# h %le\n",h);
Gp = 100.;
printf("# Gp %le\n",Gp);
epsilon = 0.01;
printf("# epsilon %le\n",epsilon);
Ep = 2.*Gp*(1.+nu);
printf("# Ep %le\n",Ep);
Em = 2.*Gp*(1.+nu)/delta;
printf("# Em %le\n",Em);

theta = gamma/(Ep*h*epsilon*epsilon);
printf("# theta %le\n",theta);
theta = 1./theta;
alpha = (Ep-Em)/(Ep+Em);
printf("# alpha %le\n",alpha);
c1 = 24.;
c2 = 18;
c3 = 6.;
c4 = -3.;
c5 = 9.;
c6 = -1.5;
for(i=0; i < 241; ++i){
j = (double) (0.01*i);
omega = 0.5*j;
denom1 =
-8.*alpha*alpha*cosh(omega)*cosh(omega)
+8.*cosh(omega)*cosh(omega)
+12.*omega*alpha
+3.*alpha*alpha*omega
-4.
+4.*alpha*alpha
+7.*sinh(omega)*cosh(omega)*alpha*alpha
+8.*cosh(omega)*sinh(omega);
tmp1 = c1*(1.-alpha)*exp(-omega)*cosh(omega)*omega/denom1;
tmp2 = c2*(1.-alpha)*exp(-omega)*alpha*cosh(omega)*cosh(omega)/denom1;
tmp3 =
-2.*cosh(omega)
-2.*sinh(omega)
-3.*cosh(omega)*alpha
+2.*sinh(omega)*alpha
+3.*omega*sinh(omega)*alpha
-3.*omega*cosh(omega)*alpha;
denom2 =
-15.*cosh(omega)*cosh(omega)*cosh(omega)*alpha*alpha
+4.*cosh(omega)
+12.*omega*cosh(omega)*alpha
+3.*cosh(omega)*alpha*alpha*omega
+11.*cosh(omega)*alpha*alpha
+15.*cosh(omega)*cosh(omega)*sinh(omega)*alpha*alpha
-12.*omega*sinh(omega)*alpha
-3.*omega*sinh(omega)*alpha*alpha
+4.*sinh(omega)
-4.*sinh(omega)*alpha*alpha;
tmp4 = c3*(1.-alpha)*exp(-omega)*cosh(omega)*tmp3/denom2;
tmp5 =
12.*omega*sinh(omega)*alpha
+3.*omega*sinh(omega)*alpha*alpha
-4.*cosh(omega)
-6.*cosh(omega)*alpha*alpha
-4.*sinh(omega)
+4.*sinh(omega)*alpha*alpha
-10.*cosh(omega)*alpha;
tmp6 = c4*cosh(omega)*tmp5/denom1;
tmp7 = c5*cosh(omega)*alpha*(4.+alpha)*omega*sinh(omega)/denom1; 
tmp8 = c6*cosh(omega)*cosh(omega)*alpha*(4.+alpha)/denom1;
f = -j*j*j + j*j*theta*(tmp1+tmp2+tmp4+tmp6+tmp7+tmp8);
printf("%le %le\n",j,f);
}
return 0;
}