/* jelolesek:   R=1          y=R-r       alpha=(Pi/2)-phi
                z=sqrt(r-rkupak)
                U(y)=V(R-y)              Uv(y)=U_yesszo(y)
                yk=ykupak                rk=rkupak
                uyk=U(yk)                uvyk=Uv(yk)
                vmin2=1-2*U(yk)
*/

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

FILE *szamok;
int uje;
long int N,i,j,n;
double alpha,eps,sinal,sin2al,dt,y,fy,uy,uvy,fvy,yuj,r,yk,
       rk,uyk,uvyk,vmin2,rk3uvyk,z,I1,I2,z2,rho,rhoinv,tmp1,tmp2,tmp3,
       IntKappa,alpha0,
       yok[101];

  /* A potencialt es a derivaltjat megado ket fuggveny */
double U(double y)
{
return((y>0)?(0.5-1/(2-log(1-y))):0);
}

double Uv(double y)
{
double r,twomlogr;
r=1-y;twomlogr=2-log(r);
return((y>0)?(1/(twomlogr*twomlogr*r)):0);
}

szamol(){

  /* nehany globalis konstans */
sinal=sin(alpha);
sin2al=sinal*sinal;
dt=sinal/N;

  /* ykupak keresese Newton-modszerrel */
y=0; fy=-sin2al; uy=0; uvy=Uv(y); fvy=2*uvy+2; yuj=y-fy/fvy; i=0;
yok[0]=y; uje=1; if(yuj==y) uje=0;
while ((i<101)&&uje) {
  y=yuj; r=1-y; uy=U(y); uvy=Uv(y);
  fy=2*r*r*uy+(2-y)*y-sin2al;
  fvy=-4*r*uy+2*r*r*uvy+2-2*y;
  yuj=y-fy/fvy;
  i++;
  yok[i]=y;
  uje=1; for(j=i-1;j&&uje;j--) if (y==yok[j]) uje=0;
  }
  if (uje||(fabs(y-yuj)>eps*fabs(y))) {
     fprintf(szamok,"Hiba: Nem konvergal a Newton-modszer.\n");
     return;
     }
yk=yuj;

  /* ertekadas a csak ykupak-tol fuggo konstansoknak */
rk=1-yk; uyk=uy; uvyk=uvy; vmin2=1-2*uyk; rk3uvyk=rk*rk*rk*uvyk;

  /* z mozgasegyenletenek numerikus megoldasa
     egyuttal ket ido szerinti integral szamitasa  */
z=sqrt(yk); I1=0; I2=0; n=0;
while(z>0){
  z2=z*z; y=yk-z2; y=((y>0)?y:0); r=rk+z2; rho=r*r; rhoinv=1/rho;
  uy=U(y); uvy=Uv(y);
  tmp1=(r+rk)*vmin2; tmp2=2*(uy-uyk);
  I1+=rhoinv;
  I2+=rhoinv*(rk3uvyk-rho*r*uvy)/(z2*tmp1-rho*tmp2);
  z-=0.5*sqrt(rhoinv*tmp1-tmp2/z2)*dt;
  n++;
  }

  /* kovetkeztetesek levonasa */
tmp3=r*uvyk/(vmin2+r*uvyk)*(1/sinal+I1*dt)+vmin2/(vmin2+r*uvyk)*I2*dt;
IntKappa-=(alpha-alpha0)*(2*sinal*tmp3-2);
fprintf(szamok,"%e\t",alpha);
fprintf(szamok,"%e\t",2*n*dt);
fprintf(szamok,"%e\t",2*I1*dt*cos(alpha));
fprintf(szamok,"%e\t",tmp3);
fprintf(szamok,"%e\n",IntKappa);
alpha0=alpha;
}

main(){

szamok=fopen("outsing.txt","w");
/* eps=1;
   while ((1+eps)!=1) eps*=0.1;
   eps*=10000; */
eps=1.0e-10;

N=100000;
fprintf(szamok,"N = %ld\n",N);

IntKappa=0; alpha0=0;

fprintf(szamok,"%11s\t%11s\t%11s\t%11s\t%11s\n","alpha","t","DeltaTheta","(2+k)/2cosfi","IntKappa");
for (alpha=0.0001;alpha<0.001;alpha+=0.0001) szamol();
for (alpha=0.001;alpha<0.01;alpha+=0.001) szamol();
for (alpha=0.01;alpha<3.141;alpha+=0.01)
szamol();

fflush(szamok);

}