#include<stdlib.h>
#include<stdio.h>
#include <string.h>
#include <stdarg.h>
#include<math.h>
#include<gsl/gsl_errno.h>
#include<gsl/gsl_odeiv.h>
#include<gsl/gsl_const_cgs.h>
#include<gsl/gsl_integration.h>
#include<gsl/gsl_matrix.h>


double numpar (char * parameter);
char * strrev(char * string);
char* stringpar (char * parameter);


/* ##################################################################################### */
/* Solution of the Poisson equation with screened potential */
/* using spherical coordinates */

/* y''[r]+2/r y[r]- k y[r]=0 */

/* The linearly independent solutions are */
/* y[r]=C1 Exp[- sqrt(k) r]/r+ C2 Exp[ sqrt(k) r]/(Sqrt[k]* r) */

/* thus there is an exponential decreasing and exponential increasing function. */
/* If one is interested only the the first one, then C2=0 */

/* Solving the equation with initial conditions requires to */
/* set the function at the value r0, y[r0] and then its derivative */
/* at the value derived analytically by derivation of the solution with C2=0 */
/* #################################################################################### */

#define Pi 3.14159



/* Jacobian */

int jac (double x, const double y[], double *dfdy, double dfdt[], void *params) {


double kpar=((double *)params)[0];
double dyi=((double *)params)[1];
double x0=((double *)params)[2];

gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy,2,2);

gsl_matrix *m=&dfdy_mat.matrix;

gsl_matrix_set(m,0,0,0);
gsl_matrix_set(m,0,1,1);
gsl_matrix_set(m,1,0,kpar);
gsl_matrix_set(m,1,1,-2/x);


dfdt[0]=0; 
dfdt[1]=2/pow(x0,2)*dyi; 

printf("Jacobian y[0] %5.3e dfdt[0] %5.3e  dfdt[1] %5.3e\n", y[0],dfdt[0], dyi);
return GSL_SUCCESS; 


}



double xmax;
char *output;

int funct (double x, const double y[], double f[], void *params); 

int main (int argc, char *argv[]) {

    FILE *dat;
  char *meth;
  double x,eps, xmin,  kpar, yi, dyi, c1, zio;

/* ######################################################################*/

if (argc < 2)  {

  printf("Test the solution of the Helmotz equation in spherical coordinates\n");
  printf("test_algorithm  yi=yi  kpar=kpar meth=meth eps=eps [xmax=xmax]\n\n");
	 exit(1);
}



/* ######################################################################*/

int t;

for (t=1; t<argc; t++) {

if (strstr(argv[t], "kpar=") != NULL)          kpar=numpar(argv[t]);
if (strstr(argv[t], "meth=") != NULL)          meth=stringpar(argv[t]);
if (strstr(argv[t], "eps=") != NULL)          eps=numpar(argv[t]);
if (strstr(argv[t], "xmax=") != NULL)          xmax=numpar(argv[t]);
if (strstr(argv[t], "yi=") != NULL)            yi=numpar(argv[t]);
/* if (strstr(argv[t], "dyi=") != NULL)          dyi=numpar(argv[t]); */
if (strstr(argv[t], "zio=") != NULL)          zio=numpar(argv[t]);

}


 
 if (!kpar) {
 printf("Parameter kpar not set\n");
 return(1);
 }

 if (! yi) {
 printf("Set the initial values for the function\n");
 return(1);
 }

if (!output) output="results_hmz.dat";
if (!eps) eps=1e-7;
if (!xmax) xmax=100.;

if (!meth) {
 printf("Parameter meth not set\n");
 return(1);
 }


 printf("kpar %5.3f  meth %s eps %5.3e\n", kpar, meth, eps);
 

/* ######################################################################*/
/* Choose the integration methd */
/* ######################################################################*/

const gsl_odeiv_step_type *T;

if (strcmp(meth,"rk4") ==0) {
/* const gsl_odeiv_step_type *T = gsl_odeiv_step_rk4; */
T = gsl_odeiv_step_rk4;


 } else if  (strcmp(meth,"gear1") ==0) {

/* const gsl_odeiv_step_type *T = gsl_odeiv_step_gear1; */
T = gsl_odeiv_step_gear1;

 } else if  (strcmp(meth,"gear2") ==0) {

/* const gsl_odeiv_step_type *T = gsl_odeiv_step_gear2; */
T = gsl_odeiv_step_gear2;


} else if  (strcmp(meth,"rkck") ==0) {

/* const gsl_odeiv_step_type *T = gsl_odeiv_step_rkck; */
T = gsl_odeiv_step_rkck;

} else if  (strcmp(meth,"rk8pd") ==0) {

/* const gsl_odeiv_step_type *T = gsl_odeiv_step_rkck; */
T = gsl_odeiv_step_rk8pd;

} else if  (strcmp(meth,"bsimp") ==0) {

/* const gsl_odeiv_step_type *T = gsl_odeiv_step_rkck; */
T = gsl_odeiv_step_bsimp;

} else {

  printf("Allowed methods are: rk4, gear1, gear2, rkck, rk8pd, bsimp\n\n"); 
  exit(1);
 }

/* void *jac; */

gsl_odeiv_step *s = gsl_odeiv_step_alloc (T,2);
gsl_odeiv_control *k = gsl_odeiv_control_y_new(eps,eps);  

gsl_odeiv_control *g = gsl_odeiv_control_yp_new(eps,eps);  
gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (2);

/* printf("Receveing r0=%5.4e  pc=%5.4e alpha=%5.4e  gamma_val=%5.4e\n", r0, pc, alpha, gamma_val); */

/* ######################################################################## */
/* Set the initial value of the independent variable and stepsize */


double xinit=pow(10,-8);
xmin=1;
x=xmin;
double h=1e-4;

/* Determine the value of the constant for y[xmin]=yi */

c1=exp(sqrt(kpar)*xmin)*xmin*yi;

/* Compute value of the derivative */

dyi=-c1* exp(-sqrt(kpar)*xmin)/pow(xmin,2)- c1* exp(-sqrt(kpar)*xmin)*sqrt(kpar)/xmin;


 dyi=-yi*(sqrt(kpar)+1/xmin);

/* ######################################################################## */

double params[3];

params[0]=kpar;
params[1]=dyi;
params[2]=xmin;

gsl_odeiv_system sys={funct, jac, 2, params};



 printf("!Value c1 %5.3f yi %5.3f  dyi %5.3e\n", c1,yi, dyi);


double y[2]={yi, dyi}; 
 
/* ######################################################################## */
/* Fixed step size integration */
/* double dydx_in[2], dydx_out[2], y_err[2]; */
/* GSL_ODEIV_FN_EVAL (&sys, t,y,dydx_in); */
/* ######################################################################## */

int status, count=0, count_loop=0;
double xprev=-10;
double minR=1000;
int countns=0;


dat=fopen(output, "w");
printf("!Value c1 %5.3f yi %5.3f  dyi %5.3e x %5.3e xmax %5.3f\n", c1,y[0], y[1], x, xmax);

while (x < xmax) {

status = gsl_odeiv_evolve_apply (e,g,s,&sys, &x,xmax , &h, y);  
 /* status = gsl_odeiv_step_apply (s,x,h,y, y_err, dydx_in, dydx_out, &sys);  */ 

 if (x > xprev) {

  fprintf(dat, "%5.4e   %5.4e\n", x, y[0]);

 }


if (status !=GSL_SUCCESS) {
break;
} 


/* ############################################################################################ */

 xprev=x;
 count++;

/* dydx_in[0]=dydx_out[0]; */
/* dydx_in[1]=dydx_out[1]; */

/*   x=x+h;  */ 
 }

 fclose(dat);

gsl_odeiv_evolve_free(e); 
gsl_odeiv_control_free(g); 
gsl_odeiv_step_free(s);

 


}


/* ######################################################################## */
/*Field equations*/
/* ######################################################################## */


int funct (double x, const double y[], double f[], void *params) {


int static iter=0, iterns=0;

double kpar=((double *)params)[0]; 

 f[0] = y[1];
 f[1] =-2 * y[1]/x + kpar* y[0];

 printf("Function x %5.3e y[0] %5.3e f[0] %5.3e\n", x, y[0], f[0]);

return GSL_SUCCESS;

 }

/* ################################################################################################ */

char * strrev(char * string) {
  int length = strlen(string);
  char * result = malloc(length+1);
  if( result != NULL ) {
    int i,j;                                         result[length] = '\0';
    for ( i = length-1, j=0;   i >= 0;   i--, j++ )  result[j] = string[i];
  }
  return result;
}

/* ################################################################################################ */

char* stringpar (char * parameter) {

int t, ch='=';
double value;
char *file;


int substr=strlen(strrev(parameter))-strlen(strchr(strrev(parameter),ch));
char *dest=malloc(substr+1);


 file=strrev((strncpy(dest, strrev(parameter), substr)));


 return(file);

}


double numpar (char * parameter) {

 

int t, ch='=';
double value;

 /*  printf("%s\n", parameter);  */


int substr=strlen(strrev(parameter))-strlen(strchr(strrev(parameter),ch));
char *dest=malloc(substr);


value=atof(strrev((strncpy(dest, strrev(parameter), substr))));


 /* printf("%s m=%lf\n", strrev(strncpy(dest, strrev(parameter), substr)), m);  */

 return(value);

}