#include #include #include #include #include #include #include #include #include #include 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 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); }