I have C code. On compliation it shows error,
warning: assignment from incompatible pointer type
Actually I wanted to evaluate integral then sum all the integrals. I would be glad if anyone could help me. Thanks
warning: assignment from incompatible pointer type
Actually I wanted to evaluate integral then sum all the integrals. I would be glad if anyone could help me. Thanks
Code:
#include <stdlib.h> #include <stdio.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_math.h> #include <gsl/gsl_vector.h> #include <gsl/gsl_multiroots.h> #include <gsl/gsl_integration.h> struct rparams { double mu; double Lambda; double G; double H; }; struct params2 { double mu; double Lambda; double H; double G; double Delta; double mue; double mu8; }; int rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f); double EQ1 (double Delta, double mue, double mu8, void *params); double EQ2 (double Delta, double mue, double mu8, void *params); double EQ3 (double Delta, double mue, double mu8, void *params); double Part1 (double p, void *params); double Part2 (double p, void *params); double Part3 (double p, void *params); double Part4 (double p, void *params); double Part5 (double p, void *params); double Part6 (double p, void *params); double Sum1 (double p, void *params); double Sum2 (double p, void *params); double Sum3 (double p, void *params); double Sum4 (double p, void *params); double Sum5 (double p, void *params); double Sum6 (double p, void *params); double function (int i, double p3, void *params, int n); double newfunction( double p3, void *params, int n); double integra_fv(int i, double Delta, double mue, double mu8, void *params); double Integrate (int i, double Delta, double mue, double mu8, void *params); double lowerlimit(double Delta, double mue, double mu8, void *params, int n); double upperlimit(double Delta, double mue, double mu8, void *params, int n); double Integralp1p2(int i, double Delta, double mue, double mu8, void *params, int n); int main (int argc, char *argv[]) { int status; double Delta, mue, mu8, mu, H; size_t iter=0; char name[1000]; const size_t n= 3; FILE *out; if(argc==3) { sprintf (name,"2SC_mu%sH%snnonzero.dat2", argv[1], argv[2]); mu=atof(argv[1]); H=atof(argv[2]); }else{ printf ("as user:\n run mu H\n"); return 1; } if ((out=fopen(name, "w"))==NULL) { printf( "Error in opening the outputfile\n"); return 1; } const gsl_multiroot_fsolver_type *T; gsl_multiroot_fsolver *s; struct rparams p= {mu, 653.0, 3.76E-6, H}; printf("H=%g mu=%g\n", H, mu); gsl_multiroot_function f= {&rosenbrock_f, n, &p}; double x_init[3]= {70.0, 150.0, 18.0}; gsl_vector *x = gsl_vector_alloc (n); gsl_vector_set (x, 0, x_init[0]); gsl_vector_set (x, 1, x_init[1]); gsl_vector_set (x, 2, x_init[2]); T= gsl_multiroot_fsolver_hybrids; s=gsl_multiroot_fsolver_alloc (T, 3); gsl_multiroot_fsolver_set (s, &f, x); // print_state (iter, s); do { iter++; printf ("iter= %3u x = % .3f % .3f % .3f" "f(x)= % .3e % .3e % .3e\n", iter, gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), gsl_vector_get (s->x, 2), gsl_vector_get (s->f, 0), gsl_vector_get (s->f, 1), gsl_vector_get (s->f, 2)); status = gsl_multiroot_fsolver_iterate (s); // print_state (iter, s); if (status) break; status = gsl_multiroot_test_residual (s->f, 1e-7); } while (status == GSL_CONTINUE && iter < 10000); printf ("status = %s\n", gsl_strerror (status)); Delta =gsl_vector_get (s->x,0); mue =gsl_vector_get (s->x, 1); mu8 =gsl_vector_get (s->x, 2); fprintf(out, "%g %g %g %g %g\n", p.mu, H, Delta, mue, mu8); fclose (out); return 0; } int rosenbrock_f (const gsl_vector * x, void *params, gsl_vector * f) { const double x0 = gsl_vector_get (x, 0); const double x1 = gsl_vector_get (x, 1); const double x2 = gsl_vector_get (x, 2); const double y0 = EQ1(x0, x1, x2, params); const double y1 = EQ2(x0, x1, x2, params); const double y2 = EQ3(x0, x1, x2, params); gsl_vector_set (f, 0, y0); gsl_vector_set (f, 1, y1); gsl_vector_set (f, 2, y2); return GSL_SUCCESS; } double EQ1 (double Delta, double mue, double mu8, void *params) { double G = ((struct rparams *) params)->G; double mu = ((struct rparams *) params)->mu; double H = ((struct rparams *) params)->H; int nmax; int n; double Pi=M_PI; double avemu= mu-1.0/6.0*mue+1.0/3.0*mu8; double P1, P2, P3, P4, f; P1= -(Delta*H)/(2.0*Pi*Pi)*Integrate(1, Delta, mue, mu8, params); /* from 0 to infty for n equal to zero */ P2= -(Delta*H)/(Pi*Pi)* integra_fv(1, Delta, mue, mu8, params); /* from 0 to infty for n not equal to zero */ P3= -(Delta*H)/(2*Pi*Pi)* Integrate(4, Delta, mue, mu8, params); /* from p2 to infty for n not equal to zero */ P4= +(Delta*H)/(2.0*Pi*Pi)*integra_fv(4, Delta, mue, mu8, params); /* from 0 to infty for n not equal to zero *****/ nmax=(avemu-sqrt(mue*mue/4.0-Delta*Delta))*(avemu-sqrt(mue*mue/4.0-Delta*Delta))/(H); double S=0; for(n=0;n<=nmax;n++){ S+=-(Delta*H)/(2.0*Pi*Pi)*Integralp1p2(1, Delta, mue, mu8, params,n)+10; } f = P1+P2+P3+P4+S+Delta/(2.0*G); printf("I'm in EQ1 nmax=%u S=%f\n", nmax, S); return f; } double EQ2 (double Delta, double mue, double mu8, void *params) { double H = ((struct rparams *) params)->H; double mu = ((struct rparams *) params)->mu; double Pi=M_PI; double mu_db, mu_ub; double P1, P2,P3,P4, f1,f; mu_db = mu +1.0/3.0*mue -2.0/3.0*mu8; mu_ub = mu-2.0/3.0*mue-2.0/3.0*mu8; P1 = H/(12.0*Pi*Pi)*Integrate(2, Delta, mue, mu8, params); /* 0 to infty for n equal to zero */ P2 = H/(6.0*Pi*Pi)*integra_fv(2, Delta, mue, mu8, params); /* 0 to infty for n not equal to zero */ P3 = H/(12.0*Pi*Pi)*integra_fv(5, Delta, mue, mu8, params); /* 0 to infty for n not equal to zero */ P4 = -H/(12.0*Pi*Pi)*Integrate(5, Delta, mue, mu8, params); /* p2 to infty for n not equal to zero */ f1 = P1 + P2+P3+P4 -(mu_db*mu_db*mu_db)/(9.0*Pi*Pi)-H/(4.0*Pi*Pi)*mue +H/(6.0*Pi*Pi)*mu_ub; int n; double s=0.0; int max = (mue*mue)/(2.0*H); for (n=1; n<=max; n++) s+=H/(2.0*Pi*Pi)*sqrt(mue*mue-2.0*n*H); int m; double f2=0.0; int maxm = (mu_ub*mu_ub)/(2.0*H); for (m=1; m<=maxm; m++) f2+=H/(3.0*Pi*Pi)*sqrt(mu_ub*mu_ub-2.0*m*H); f=f1+f2-s; /* printf("I am in eq2 f2=%f mue=%f mu_db=%f mu_ub=%f maxm=%u\n", f2, mue, mu_db, mu_ub, maxm); */ return f; } double EQ3(double Delta, double mue, double mu8, void *params) { double H = ((struct rparams *) params)->H; double mu = ((struct rparams *) params)->mu; double Pi=M_PI; double mu_db = mu +1.0/3.0*mue -2.0/3.0*mu8; double mu_ub = mu-2.0/3.0*mue-2.0/3.0*mu8; double P1, P2, P3, P4, f; P1 = -H/(6.0*Pi*Pi)*Integrate(3, Delta, mue, mu8, params); /* 0 to infty for n equal to zero */ P2 = -H/(3.0*Pi*Pi)*integra_fv(3, Delta, mue, mu8, params); /* 0 to infty for n not equal to zero */ P3 = -H/(6.0*Pi*Pi)*integra_fv(6, Delta, mue, mu8, params); /* 0 to infty for n not equal to zero */ P4 = H/(6.0*Pi*Pi)*Integrate(6, Delta, mue, mu8, params); /* p2 to infty for n not equal to zero */ int j; double f3=0; int maxmm =mu_ub*mu_ub/(2.0*H); for(j=1; j<=maxmm; j++) f3+=H/(3.0*Pi*Pi)*sqrt(mu_ub*mu_ub-2.0*j*H); f = P1 + P2 +P3+P4+ f3 + 2.0*(mu_db*mu_db*mu_db)/(9.0*Pi*Pi)+H/(6.0*Pi*Pi)*mu_ub; return f; } double Part1 (double p3, void *params) { double Delta = ((struct params2 *) params)->Delta; double mu = ((struct params2 *) params)->mu; double mue = ((struct params2 *) params)->mue; double mu8 = ((struct params2 *) params)->mu8; double avemu; double P1; avemu = mu-mue/6.0+mu8/3.0; P1 = (1.0/sqrt((p3+avemu)*(p3+avemu)+Delta*Delta)+1.0/(sqrt((p3-avemu)*(p3-avemu)+Delta*Delta))); return P1; } /* first part for mue equation*/ double Part2 (double p3, void *params) { double Delta = ((struct params2 *) params)->Delta; double mu = ((struct params2 *) params)->mu; double mue = ((struct params2 *) params)->mue; double mu8 = ((struct params2 *) params)->mu8; double avemu; double P2; avemu = mu-mue/6.0+mu8/3.0; P2 = ((p3+avemu)/(sqrt((p3+avemu)*(p3+avemu)+Delta*Delta))-(p3-avemu)/(sqrt((p3-avemu)*(p3-avemu)+Delta*Delta))); return P2; } /* first part for mu8 equation*/ double Part3 (double p3, void *params) { double Delta = ((struct params2 *) params)->Delta; double mu = ((struct params2 *) params)->mu; double mue = ((struct params2 *) params)->mue; double mu8 = ((struct params2 *) params)->mu8; double avemu; double P3; avemu = mu-mue/6.0+mu8/3.0; P3 = ((p3+avemu)/sqrt((p3+avemu)*(p3+avemu)+Delta*Delta)-(p3-avemu)/sqrt((p3-avemu)*(p3-avemu)+Delta*Delta)); /*printf("P3=%f\n", P3); */ return P3; } /* part containing integration from p2 to infty */ double Part4 (double p3, void *params) { double Delta = ((struct params2 *) params)->Delta; double Lambda = ((struct params2 *) params)->Lambda; double mu = ((struct params2 *) params)->mu; double mue = ((struct params2 *) params)->mue; double mu8 = ((struct params2 *) params)->mu8; double P4; double avemu = mu-mue/6.0+mu8/3.0; double b1=sqrt(p3*p3+avemu*avemu); if(mue/2.0>Delta){ P4=1.0/sqrt((b1-avemu)*(b1-avemu)+Delta*Delta); /* second part mue/2 greater than Delta for Delta equation */ }else{ P4=0; } return P4*exp(-avemu*avemu/(Lambda*Lambda)); } double Part5 (double p3, void *params) { double Delta = ((struct params2 *) params)->Delta; double Lambda = ((struct params2 *) params)->Lambda; double mu = ((struct params2 *) params)->mu; double mue = ((struct params2 *) params)->mue; double mu8 = ((struct params2 *) params)->mu8; double P5; double avemu = mu-mue/6.0+mu8/3.0; double b1=sqrt(p3*p3+avemu*avemu); if(mue/2.0>Delta){ P5=((b1-avemu)/sqrt((b1-avemu)*(b1-avemu)+Delta*Delta)-3.0); /* second part mue/2 greater than Delta for mue equation */ }else{ P5=0; } return P5*exp(-avemu*avemu/(Lambda*Lambda)); } double Part6 (double p3, void *params) { double Delta = ((struct params2 *) params)->Delta; double Lambda = ((struct params2 *) params)->Lambda; double mu = ((struct params2 *) params)->mu; double mue = ((struct params2 *) params)->mue; double mu8 = ((struct params2 *) params)->mu8; double P6; double avemu = mu-mue/6.0+mu8/3.0; double b1=sqrt(p3*p3+avemu*avemu); if(mue/2.0>Delta){ P6=(b1-avemu)/sqrt((b1-avemu)*(b1-avemu)+Delta*Delta); /* second part mue/2 greater than Delta for mue equation */ }else{ P6=0; } return P6*exp(-avemu*avemu/(Lambda*Lambda)); } double Sum1(double p, void *params) { double Lambda = ((struct params2 *) params)->Lambda; double H = ((struct params2 *) params)->H; double S=0; int n; int max = (10.0*Lambda*Lambda)/H; for (n=1;n<=max; n++){ S+=function(1,p, params, n); } return S; } double Sum2(double p, void *params) { double Lambda = ((struct params2 *) params)->Lambda; double H = ((struct params2 *) params)->H; double S=0; int n; int max = (10.0*Lambda*Lambda)/H; for (n=1;n<=max; n++){ S+=function(2,p, params, n); } return S; } double Sum3(double p, void *params) { double Lambda = ((struct params2 *) params)->Lambda; double H = ((struct params2 *) params)->H; double S=0; int n; int max = (10.0*Lambda*Lambda)/H; for (n=1;n<=max; n++){ S+=function(3,p, params, n); } return S; } /* for mue/2 greater than Delta * * * */ double Sum4(double p, void *params) { double mu =((struct params2*) params)->mu; double H = ((struct params2*) params)->H; double Lambda = ((struct params2*) params)->Lambda; double mue =((struct params2*) params)->mue; double mu8 =((struct params2*) params)->mu8; double avemu = mu-mue/6.0+mu8/3.0; double S=0; int n; int max = 10.0*Lambda*Lambda/H; for (n=avemu*avemu/H;n<=max; n++){ S+=function(4,p, params, n); } return S; } double Sum5(double p, void *params) { double mu =((struct params2*) params)->mu; double H = ((struct params2*) params)->H; double Lambda = ((struct params2*) params)->Lambda; double mue =((struct params2*) params)->mue; double mu8 =((struct params2*) params)->mu8; double avemu = mu-mue/6.0+mu8/3.0; double S=0; int n; int max = 10.0*Lambda*Lambda/H; for (n=avemu*avemu/H;n<=max; n++){ S+=function(5,p, params, n); } return S; } double Sum6(double p, void *params) { double mu =((struct params2*) params)->mu; double H = ((struct params2*) params)->H; double Lambda = ((struct params2*) params)->Lambda; double mue =((struct params2*) params)->mue; double mu8 =((struct params2*) params)->mu8; double avemu = mu-mue/6.0+mu8/3.0; double S=0; int n; int max = 10.0*Lambda*Lambda/H; for (n=avemu*avemu/H;n<=max; n++){ S+=function(6,p, params, n); } return S; } double function(int i, double p3, void *params, int n) { double mu = ((struct params2 *) params )->mu; double Lambda = ((struct params2 *) params)->Lambda; double Delta = ((struct params2 *) params)->Delta; double mue = ((struct params2 *) params)->mue; double mu8 = ((struct params2 *) params)->mu8; double H = ((struct params2 *) params)->H; double a1, avemu, f1, f2, nint; avemu = mu-mue/6.0+mu8/3.0; nint =n*1.0; a1=sqrt(p3*p3+nint*H); if(i==1){ f1=1.0/sqrt((a1+avemu)*(a1+avemu)+Delta*Delta); /* first part for Delta equation */ f2=1.0/sqrt((a1-avemu)*(a1-avemu)+Delta*Delta); return (f1+f2)*exp((-nint*H)/(Lambda*Lambda)); }else{ if(i==2){ f1=(a1+avemu)/sqrt((a1+avemu)*(a1+avemu)+Delta*Delta); /* first part for mue equation */ f2=(a1-avemu)/sqrt((a1-avemu)*(a1-avemu)+Delta*Delta); return (f1-f2)*exp((-nint*H)/(Lambda*Lambda)); }else{ if(i==3){ f1=(a1+avemu)/sqrt((a1+avemu)*(a1+avemu)+Delta*Delta); /* first part of mu8 equation */ f2=(a1-avemu)/sqrt((a1-avemu)*(a1-avemu)+Delta*Delta); return (f1-f2)*exp((-nint*H)/(Lambda*Lambda)); }else{ if(i==4){ if(mue/2.0>=Delta){ f1=1.0/sqrt((a1-avemu)*(a1-avemu)+Delta*Delta); /* second part mue/2 greater than Delta for Delta equation */ }else{ f1=0; } return f1*exp((-nint*H)/(Lambda*Lambda)); }else{ if(i==5){ if(mue/2.0>=Delta){ f1=((a1-avemu)/sqrt((a1-avemu)*(a1-avemu)+Delta*Delta)-3.0); /* second part mue/2 greater than Delta for mue equation */ }else{ f1=0; } return f1*exp((-nint*H)/(Lambda*Lambda)); }else{ if(mue/2.0>=Delta){ f1=(a1-avemu)/sqrt((a1-avemu)*(a1-avemu)+Delta*Delta); /* second part mue/2 greater than Delta for mu8 equation */ }else{ f1=0; } return f1*exp((-nint*H)/(Lambda*Lambda)); } } } } } } double newfunction( double p3, void *params, int n) { double mu = ((struct params2 *) params )->mu; double Lambda = ((struct params2 *) params)->Lambda; double Delta = ((struct params2 *) params)->Delta; double mue = ((struct params2 *) params)->mue; double mu8 = ((struct params2 *) params)->mu8; double H = ((struct params2 *) params)->H; double a1, avemu, f1, nint; avemu = mu-mue/6.0+mu8/3.0; nint =n*1.0; a1=sqrt(p3*p3+nint*H); if(mue/2.0>=Delta){ f1=1.0/sqrt((a1-avemu)*(a1-avemu)+Delta*Delta); /* second part mue/2 greater than Delta for Delta equation */ }else{ f1=0; } return f1*exp((-nint*H)/(Lambda*Lambda)); } double integra_fv(int i, double Delta, double mue, double mu8, void *params) /* integration for n not equal to zero */ { gsl_integration_workspace * w = gsl_integration_workspace_alloc (10000); double result, error; double mu = ((struct rparams *) params)->mu; double Lambda = ((struct rparams *) params)->Lambda; double H = ((struct rparams *) params)->H; double G= ((struct rparams *) params)->G; struct params2 p = {mu, Lambda, H, G, Delta, mue, mu8}; gsl_function F; if(i==1){ F.function = &Sum1; F.params = &p; }else{ if(i==2){ F.function = &Sum2; F.params = &p; }else{ if(i==3){ F.function = &Sum3; F.params = &p; }else{ if(i==4){ F.function = &Sum4; F.params = &p; }else{ if(i==5){ F.function = &Sum5; F.params = &p; }else{ F.function = &Sum6; F.params = &p; } } } } } gsl_integration_qag(&F, 0.0, 653.0, 0.0, 1e-7, 10000, 6, w, &result, &error); gsl_integration_workspace_free(w); return result; } /* integration for n equal to zero */ double Integrate(int i, double Delta, double mue, double mu8, void *params) { gsl_integration_workspace * w = gsl_integration_workspace_alloc(10000); double result, error; double mu = ((struct rparams *) params)->mu; double Lambda = ((struct rparams *) params)->Lambda; double H = ((struct rparams *) params)->H; double G= ((struct rparams *) params)->G; struct params2 p = {mu, Lambda, H, G, Delta, mue, mu8}; gsl_function F; if(i==1){ F.function = &Part1; F.params = &p; }else{ if(i==2){ F.function = &Part2; F.params = &p; }else{ if(i==3){ F.function = &Part3; F.params = &p; }else{ if(i==4){ F.function = &Part4; F.params = &p; }else{ if(i==5){ F.function = &Part5; F.params = &p; }else{ F.function = &Part6; F.params = &p; } } } } } gsl_integration_qag (&F, 0.0, 653.0, 0.0, 1e-7, 10000, 6, w, &result, &error); gsl_integration_workspace_free (w); return result; } double lowerlimit(double Delta, double mue, double mu8, void *params, int n) { double mu = ((struct params2 *) params)->mu; double H = ((struct params2 *) params)->H; double S; double nint=1.0*n; double avemu=mu-mue*1.0/6.0+mu8*1.0/3.0; if(mue/2.0>Delta){ S=sqrt((avemu-sqrt(mue*mue/4.0-Delta*Delta))*(avemu-sqrt(mue*mue/4.0-Delta*Delta))-nint*H); }else{ S=0; } /* printf("the variables are mue=%f avemu=%f mu=%f s=%f\n", mue, avemu, mu, S); */ /* printf("lowerlimit=%f\n", S); */ return S; } double upperlimit(double Delta, double mue, double mu8, void *params, int n) { double mu = ((struct params2 *) params)->mu; double H = ((struct params2 *) params)->H; double S; double nint=1.0*n; double avemu=mu-mue/6.0+mu8*1.0/3.0; if(mue/2.0>Delta){ S=sqrt((avemu+sqrt(mue*mue/4.0-Delta*Delta))*(avemu+sqrt(mue*mue/4.0-Delta*Delta))-nint*H); }else{ S=0; } /* printf("upperlimit S=%f\n", S); */ return S; } double Integralp1p2(int i, double Delta, double mue, double mu8, void *params, int n) { gsl_integration_workspace * w = gsl_integration_workspace_alloc (10000); double result, error; double mu = ((struct rparams *) params)->mu; double Lambda = ((struct rparams *) params)->Lambda; double H = ((struct rparams *) params)->H; double G= ((struct rparams *) params)->G; struct params2 p = {mu, Lambda, H, G, Delta, mue, mu8}; gsl_function F; if(i==1){ F.function = &newfunction; F.params = &p; }else{ if(i==2){ F.function = &newfunction; F.params = &p; }else{ if(i==3){ F.function = &newfunction;; F.params = &p; }else{ F.function = &newfunction;; F.params = &p; } } } double p1=lowerlimit(Delta, mue, mu8, params,n); double p2=upperlimit(Delta, mue, mu8, params,n); gsl_integration_qag(&F, p1, p2, 0.0, 1e-7, 10000, 6, w, &result, &error); gsl_integration_workspace_free(w); return result; }
Comment