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

#define SIZE 2
#define wmax 2
#define wmin 0.001
#define rhor 2.0001
#define rinf 1000.
#define eps 1


double r;
double complex X[SIZE],X1[SIZE];
double rh,ri,h,wi,wf,e,hw;
int    finesse,i,j,k,kl,l,n;
// DE system to integrate (n=2)
      void F(double r, double w, int l, double f,double e, double complex x[], double complex xp[]) {
        xp[0] = x[1];
        xp[1] = ((-(2-r)*(2*r*r+l*(l+1)*r*r*r + l*(l+1)*e) - r*r*r*r*r*r*w*w)*x[0]+2*(2-r)*r*r*r*x[1])/(r*r*r*r*(r-2)*(r-2));
      }
//Integrate first order DE sytem from x to x+h using Runge-Kutta
      void RK4n(int n, double r, double w, int l, double f,double e, double complex x[], double h, double complex x1[]) {
        double complex c1[SIZE],c2[SIZE],c3[SIZE],c4[SIZE],xx[SIZE];
        double h2;
        int i;
        F(r,w,l,f,e,x,c1);
        h2=h/2.0;
        for (i=0; i<n; i++) xx[i]=x[i]+h2*c1[i];
        F(r+h2,w,l,f,e,xx,c2);
        for (i=0; i<n; i++) xx[i]=x[i]+h2*c2[i];
        F(r+h2,w,l,f,e,xx,c3);
        for (i=0; i<n; i++) xx[i]=x[i]+h*c3[i];
        F(r+h,w,l,f,e,xx,c4);
        for (i=0; i<n; i++)
          x1[i]=x[i]+h*(c1[i]+2.0*c2[i]+2.0*c3[i]+c4[i])/6.0;
      }
int main(){

double w,f;
double complex dXi;
n=2;           // size of DE system
e=eps;          //value of epsilon
rh = rhor;      //starting r
ri = rinf;      //ending r
wi=wmin;        // starting w
wf=wmax;        // ending w
kl=25;         // number of calculated points
finesse=20;
h=0.1;
r=rh;
w=wi;
f=(1-2/r);         //f(r) function of the metric
dXi=(-I*w)/(1-2/rh);
double complex tr2;
l=0;

X[0]=1;    // starting X values
X[1]=dXi;  // starting X' values

FILE *fpout;
fpout = fopen("Transe0l0c.dat","w");
// integration loop
while (w<wmax){
    while(r<ri){
        RK4n(n,r,w,l,f,e,X,h,X1);
        for (j=0; j<n; j++) X[i]=X1[i];
        r+=0.1;
        }
    tr2=4/(2+X[0]*conj(X[0])+((1-2/ri)*(1-2/ri)/(w*w))*X[1]*conj(X[1]));
    fprintf(fpout, "%.8lf %.8lf\n",w,creal(tr2));
    w+=0.01;
    X[0]=1;
    X[1]=dXi;
    }

fclose(fpout);

return (0);
}