source: Sophya/trunk/Cosmo/SimLSS/hu_sigma8.c@ 3495

Last change on this file since 3495 was 3315, checked in by cmv, 18 years ago

add code fct transfert de W.Hu + prog de comparaison avec pkspectrum.cc cmv 24/08/2007
File size: 18.7 KB
Line 
1/* Code to calculate sigma8 for COBE-normalized CDM and MDM models,
2using the approximation of Eisenstein & Hu (1997). */
3
4/* This code will calculate sigma_8 and other quantities for CDM and
5MDM models. Because the fitting form for the transfer function does not
6include the baryon oscillations, the answer will not be accurate for
7large baryon fractions. The COBE normalization is from the papers
8on Bunn & White (1997) and Hu & White (1997). */
9
10/* Accuracy is generally better than 5%. Recall that the COBE normalization
11is only 7% accurate (1-sigma). This program gives low answers compared to
12CMBfast for lower values of Omega_0*h^2, and vice versa for higher values.
13Values of Omega_0 near 0.2 can produce answers that are 5% low. */
14
15/* The integrator used here is a bit crude, but suffices for rapid
1610^-3 accuracy. You could get better performance by substituting an
17integrator from a numerical package. See the end of the file for details. */
18
19/* This program uses the routines TFmdm_set_cosm() and TFmdm_onek_hmpc(),
20available from http://www.sns.ias.edu/~whu/power/power.html in the file
21power.c. To compile this program, you'll need that file as well, e.g.
22 cc -O sigma8.c power.c -o sigma8 -lm */
23
24/* Lists of versions and bugs:
25v1.0 -- 10/23/97 -- Original release version.
26v1.1 -- 11/13/97 -- We correct a missing factor of omega_m in the
27 conversion between mass and length scales [lines
28 32 and 34 of main()]. This only affects the results
29 if you explicitly used the -mass option.
30v1.2 -- 04/08/98 -- We correct a typo in the second-order tilt scaling
31 of the COBE normalization for flat models with tensors.
32v1.3 -- 09/15/98 -- We correct a typo in the tilt scaling of the COBE
33 normalization for open models with tensors.
34v1.4 -- 12/02/98 -- We correct the COBE normalization for tilted open models,
35 undoing the change of v1.3 and fixing what was
36 wrong there to begin with.
37*/
38
39#include <math.h>
40#include <stdio.h>
41#include <stdlib.h>
42#include <string.h>
43
44int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm,
45 int degen_hdm, double omega_lambda, double hubble, double redshift);
46double TFmdm_onek_hmpc(double kk);
47
48double sig2tophat(double x);
49double sig2gauss(double x);
50double sig2sharpk(double x);
51double sig2tophatnu(double x);
52double sig2gaussnu(double x);
53double sig2sharpknu(double x);
54
55#define WEIRD -2538 /* A strange number, unlikely to be entered */
56#define RHOCRIT 2.78e11 /* The critical density, for h=1, in Msun/Mpc^3 */
57
58/* Convenience from Numerical Recipes in C, 2nd edition
59------------ COMPLETEMENT DEBILE (cmv) ????????? ------------
60static double sqrarg;
61#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
62*/
63#define SQR(a) (a*a)
64
65/* Global variable */
66double omega_m, omega_b, omega_n,
67 numdegen, lambda, redshift,
68 hubble, scale, mass, tilt;
69double deltaH; /* COBE-normalization */
70double sig2_scale; /* For passing a length scale to the integrator */
71int qtensors, quiet, qredshift;
72
73/* ========================== User interaction routines ============== */
74void usage()
75/* Print out information about the allowed options and then quit. */
76{
77 fprintf(stderr,"Command line options:\n");
78 fprintf(stderr,"-om or -omega <omega_matter, including CDM, baryons, massive neutrinos>\n");
79 fprintf(stderr,"-ob <omega_baryon> -fb <baryon_fraction>\n");
80 fprintf(stderr,"-on <omega_neutrino> -fn <neutrino_fraction>\n");
81 fprintf(stderr,"-numdegen or -Nnu <number of degenerate neutrino species>, default 1\n");
82 fprintf(stderr,"-lambda <omega_lambda> -flat (choose lambda appropriately)\n");
83 fprintf(stderr,"-h or -hubble <H_0, units of 100 km/s/Mpc>\n");
84 fprintf(stderr,"-H or -Hubble <H_0, units of 1 km/s/Mpc>\n");
85 fprintf(stderr,"-tilt or -n <scalar spectral tilt>, no tensors.\n");
86 fprintf(stderr,"-Tilt or -N <scalar spectral tilt>, power-law tensors.\n");
87 fprintf(stderr,"-z <redshift>, if none given then appropriate value assumed.\n");
88 fprintf(stderr,"The following request the rms fluctuations inside a tophat window function:\n");
89 fprintf(stderr,"-scale <comoving length scale, h^-1 Mpc>\n");
90 fprintf(stderr,"-mass <mass scale, h^-1 Msun>\n");
91 fprintf(stderr,"If neither -scale or -mass is given, sigma_8 and other quantities are computed.\n");
92 fprintf(stderr,"-quiet (print only answer to stdout)\n");
93 exit(1);
94 return;
95}
96
97void crash(char *s)
98{
99 fprintf(stderr,"%s\n", s); exit(1); return;
100}
101
102void warn(char *s)
103{
104 fprintf(stderr,"%s\n", s); return;
105}
106
107
108void parse_command_line(int argc, char *argv[])
109/* Parse the command line, and look for illegal entries */
110{
111 int i;
112 if (argc==1) usage();
113 i = 1;
114 omega_m = 1.0;
115 omega_b = 0.05;
116 omega_n = 0.0;
117 numdegen = 1.0;
118 lambda = 0.0;
119 redshift = 0.0;
120 hubble = 0.5;
121 scale = -1.0;
122 mass = -1.0;
123 qtensors = 0;
124 quiet = 0;
125 tilt = 1.0;
126 qredshift = 0;
127
128 while (i<argc) {
129 if (!strcmp(argv[i],"-om")||!strcmp(argv[i],"-omega")) {
130 ++i;
131 if (i>=argc) usage();
132 sscanf(argv[i],"%lf",&omega_m); ++i;
133 } else if (!strcmp(argv[i],"-ob")) {
134 ++i;
135 if (i>=argc) usage();
136 sscanf(argv[i],"%lf",&omega_b); ++i;
137 } else if (!strcmp(argv[i],"-fb")) {
138 ++i;
139 if (i>=argc) usage();
140 sscanf(argv[i],"%lf",&omega_b); ++i;
141 omega_b *= -1.0; /* Signal that a fraction was desired */
142 } else if (!strcmp(argv[i],"-on")) {
143 ++i;
144 if (i>=argc) usage();
145 sscanf(argv[i],"%lf",&omega_n); ++i;
146 } else if (!strcmp(argv[i],"-fn")) {
147 ++i;
148 if (i>=argc) usage();
149 sscanf(argv[i],"%lf",&omega_n); ++i;
150 omega_n *= -1.0; /* Signal that a fraction was desired */
151 } else if (!strcmp(argv[i],"-h")||!strcmp(argv[i],"-hubble")) {
152 ++i;
153 if (i>=argc) usage();
154 sscanf(argv[i],"%lf",&hubble); ++i;
155 } else if (!strcmp(argv[i],"-H")||!strcmp(argv[i],"-Hubble")) {
156 ++i;
157 if (i>=argc) usage();
158 sscanf(argv[i],"%lf",&hubble); ++i;
159 hubble /= 100.0; /* Units were 1 km/s/Mpc */
160 } else if (!strcmp(argv[i],"-lambda")) {
161 ++i;
162 if (i>=argc) usage();
163 sscanf(argv[i],"%lf",&lambda); ++i;
164 } else if (!strcmp(argv[i],"-flat")) {
165 ++i;
166 lambda = WEIRD;
167 } else if (!strcmp(argv[i],"-numdegen")||!strcmp(argv[i],"-Nnu")) {
168 ++i;
169 if (i>=argc) usage();
170 sscanf(argv[i],"%lf",&numdegen); ++i;
171 } else if (!strcmp(argv[i],"-n")||!strcmp(argv[i],"-tilt")) {
172 ++i;
173 if (i>=argc) usage();
174 sscanf(argv[i],"%lf",&tilt); ++i;
175 qtensors = 0;
176 } else if (!strcmp(argv[i],"-N")||!strcmp(argv[i],"-Tilt")) {
177 ++i;
178 if (i>=argc) usage();
179 sscanf(argv[i],"%lf",&tilt); ++i;
180 qtensors = 1;
181 } else if (!strcmp(argv[i],"-z")) {
182 ++i;
183 if (i>=argc) usage();
184 sscanf(argv[i],"%lf",&redshift); ++i;
185 qredshift = 1;
186 } else if (!strcmp(argv[i],"-scale")) {
187 ++i;
188 if (i>=argc) usage();
189 sscanf(argv[i],"%lf",&scale); ++i;
190 } else if (!strcmp(argv[i],"-mass")) {
191 ++i;
192 if (i>=argc) usage();
193 sscanf(argv[i],"%lf",&mass); ++i;
194 } else if (!strcmp(argv[i],"-quiet")) {
195 ++i;
196 quiet = 1;
197 } else usage();
198 }
199
200 /* Look for illegal entries */
201 if (omega_m<=0.0) crash("Omega must be positive.");
202 if (hubble<=0.0) crash("Hubble constant must be positive.");
203 else if (hubble>10.0) warn("Did you get the Hubble units right?");
204 if (numdegen<0) {
205 warn("Negative number of degenerate species reset to 1.");
206 numdegen = 1.0;
207 }
208 if (redshift<=-1.0) crash("Redshift must exceed -1.");
209 if (scale>0.0 && mass>0.0) {
210 warn("Input length OR mass scale. Defaulting to length.");
211 mass = -1.0;
212 }
213
214 /* Now look for fractions, if they were input. */
215 if (omega_b<0.0) omega_b = omega_m*(-omega_b);
216 if (omega_n<0.0) omega_n = omega_m*(-omega_n);
217 if (lambda== WEIRD) lambda = 1-omega_m; /* Set to be flat */
218
219 if (omega_b+omega_n>omega_m) crash("Baryons + Neutrinos exceeds Omega_0.");
220 else if (omega_b+omega_n>0.6*omega_m)
221 warn("CDM fraction is below range of TF fitting formula.");
222 if (tilt>1.3 || tilt < 0.7)
223 warn("Tilt out of range of Bunn and White fitting formula.");
224 if (tilt!=1 && (omega_m<0.2 || omega_m>1.0))
225 warn("Omega out of range of Bunn and White fitting formula.");
226 if (tilt==1 && (omega_m<0.2 || omega_m>1.6 || lambda<0.0 || lambda>0.8))
227 warn("Omega or Lambda out of range of Bunn and White fitting formula.");
228 if (qtensors!=0&&(tilt>1.0||
229 (fabs(lambda)>1e-5&&fabs(1.0-omega_m-lambda)>1e-5)))
230 crash("Tensor formulae only hold for tilts < 1 and simple cosmologies.");
231 if (omega_b/omega_m>0.4) warn("Baryon oscillations will be large. TF fitting form will be inaccurate.");
232 if (omega_n/omega_m>0.4) warn("Neutrino fraction outside of range of TF fitting form.");
233 if (omega_m*hubble*hubble>0.4) warn("Omega_0 h^2 outside of range of TF fitting form.");
234 if (omega_m*hubble*hubble<0.06) warn("Omega_0 h^2 outside of range of TF fitting form.");
235 if (redshift>29) warn("Redshift is higher than range of TF fitting form.");
236
237 return;
238}
239
240/* ==================================== MAIN() ========================= */
241
242void main(int argc, char *argv[])
243{
244 double erfcc(double x);
245 double sigma, sigma8;
246 extern double growth_to_z0, z_equality, alpha_nu, beta_c;
247 double cobenorm();
248 double integrate(double (*func)(double));
249 parse_command_line(argc, argv);
250
251 /* Set all the background quantities */
252 TFmdm_set_cosm(omega_m, omega_b, omega_n, (int) numdegen, lambda,
253 hubble, redshift);
254 deltaH = cobenorm();
255
256 if (!quiet) { /* Echo cosmological information to the screen */
257 printf("Omega_0 = %5.3f, Omega_baryon = %5.3f, Omega_neutrino = %5.3f\n",
258 omega_m, omega_b, omega_n);
259 printf("Hubble = %5.3f, Lambda = %5.3f, N_neutrinos = %d\n",
260 hubble, lambda, (int) numdegen);
261 printf("Redshift = %5.3f, Growth = %6.4f relative to z=0\n",
262 redshift, growth_to_z0);
263 printf("Tilt = %5.3f, COBE normalization: deltaH = %8.3e",
264 tilt, deltaH);
265 if (tilt!=1.0)
266 if (qtensors) printf(" with tensors.\n");
267 else printf(" without tensors.\n");
268 else printf("\n");
269 }
270 /* If we've been asked to compute for a particular scale, do that */
271 if (scale>0 || mass>0) {
272 if (mass<0)
273 mass = 4.0*M_PI/3.0*omega_m*RHOCRIT*scale*scale*scale;
274 else if (scale<0)
275 scale = pow(3.0/4.0/M_PI/RHOCRIT/omega_m*mass, 1.0/3.0);
276 /* Integrate tophat for the given scale */
277 sig2_scale = scale;
278 sigma = sqrt(integrate(sig2tophatnu));
279
280 /* Print out the result */
281 if (!quiet) printf("For length scale %6.4f h^-1 Mpc, mass scale %8.3e h^-1 Msun.\nrms tophat fluctuations are ",
282 sig2_scale, mass);
283 printf("%6.4f\n", sigma);
284 } else {
285 /* No scale was requested, so we do the standard tasks */
286 /* Calculate sigma_8 */
287 sig2_scale = 8.0;
288 sigma8 = sqrt(integrate(sig2tophatnu));
289 if (quiet) {
290 printf("%6.4f\n", sigma8);
291 return; /* This is all we're doing */
292 }
293 printf("Sigma_8 = %6.4f\n", sigma8);
294
295 /* Calculate sigma_50 */
296 sig2_scale = 50.0;
297 sigma = sqrt(integrate(sig2tophatnu));
298 printf("Sigma_50 = %6.4f Sigma_50/sigma_8 = %6.4f\n",
299 sigma, sigma/sigma8);
300
301 /* We default to z=4 for DLA */
302 if (!qredshift) {
303 redshift = 4.0;
304 TFmdm_set_cosm(omega_m, omega_b, omega_n, (int) numdegen, lambda,
305 hubble, redshift);
306 }
307 sig2_scale = 50.0/100.0/sqrt(1+redshift)/sqrt(2*3.14159)*
308 pow(8.0*3.14159/3.0/sqrt(89.0),1./3.);
309 sigma = sqrt(integrate(sig2gaussnu));
310 printf("Omega_gas = %8.3e at redshift %5.3f due to \n sigma(50 km/s) = sigma(%5.3f h^-1 Mpc) = %6.4f (DLA, Klypin et al. 1994)\n",
311 omega_b*erfcc(1.33/1.414/sigma), redshift, sig2_scale, sigma);
312
313 /* We default to z=3 for Lyman alpha forest */
314 if (!qredshift) {
315 redshift = 3.0;
316 TFmdm_set_cosm(omega_m, omega_b, omega_n, (int) numdegen, lambda,
317 hubble, redshift);
318 }
319 sig2_scale = 1.0/(34.0*sqrt(omega_m/2.0)); /* Gnedin's k_34 */
320 sigma = sqrt(integrate(sig2gauss));
321 printf("Sigma_Lya = %6.4f at redshift %5.3f (Ly-alpha forest, Gnedin 1997)\n",
322 sigma, redshift);
323
324 printf("Small-scale: alpha_nu = %6.4f, beta_c = %6.4f\n",
325 alpha_nu, beta_c);
326 }
327 return;
328}
329
330/* ========================== Integration Kernals ======================= */
331
332double W2tophat(double x)
333/* The square of the tophat window */
334{
335 double j1onx;
336 if (x<0.03)
337 j1onx = (1.0/3.0-x*x*1.0/30.0);
338 else j1onx = (cos(x)-sin(x)/x)/x/x;
339 return 9.0*j1onx*j1onx;
340}
341
342double sig2tophat(double x)
343/* This is the integrand for the mass fluctuation integral
344sig2 = \int_0^\infty (dk/k) Delta^2(k) W^2(kr)
345broken and recombined so that one integrates on (0,1] in a new variable x. */
346/* Here, the window is a real-space tophat of radius sig2_scale, the
347latter being passed from outside in units of h^-1 Mpc */
348{
349 double powerspcb(double k);
350 return (powerspcb(x/sig2_scale)*W2tophat(x)+
351 powerspcb(1.0/x/sig2_scale)*W2tophat(1.0/x))/x;
352}
353
354double sig2tophatnu(double x)
355/* Same as sig2tophat(), but with P_cbnu */
356{
357 double powerspcbnu(double k);
358 return (powerspcbnu(x/sig2_scale)*W2tophat(x)+
359 powerspcbnu(1.0/x/sig2_scale)*W2tophat(1.0/x))/x;
360}
361
362double sig2gauss(double x)
363/* This is the integrand for the mass fluctuation integral
364sig2 = \int_0^\infty (dk/k) Delta^2(k) W^2(kr)
365broken and recombined so that one integrates on (0,1] in a new variable x. */
366/* Here, the window is a gaussian of radius sig2_scale, the
367latter being passed from outside in units of h^-1 Mpc */
368{
369 double powerspcb(double k);
370 return (powerspcb(x/sig2_scale)*exp(-x*x)+
371 powerspcb(1.0/x/sig2_scale)*exp(-1.0/x/x))/x;
372}
373
374double sig2gaussnu(double x)
375/* Same as sig2gauss(), but with P_cbnu */
376{
377 double powerspcbnu(double k);
378 return (powerspcbnu(x/sig2_scale)*exp(-x*x)+
379 powerspcbnu(1.0/x/sig2_scale)*exp(-1.0/x/x))/x;
380}
381
382double sig2sharpk(double x)
383/* This is the integrand for the mass fluctuation integral
384sig2 = \int_0^\infty (dk/k) Delta^2(k) W^2(kr)
385broken and recombined so that one integrates on (0,1] in a new variable x. */
386/* Here, the window is a k-space tophat of radius 1/sig2_scale, the
387latter being passed from outside in units of h^-1 Mpc */
388{
389 double powerspcb(double k);
390 return (powerspcb(x/sig2_scale))/x;
391}
392
393double sig2sharpknu(double x)
394/* Same as sig2sharpk(), but with P_cbnu */
395{
396 double powerspcbnu(double k);
397 return (powerspcbnu(x/sig2_scale))/x;
398}
399
400/* ============================ Power Spectrum and COBE ================ */
401double powerspcbnu(double k)
402/* Returns Delta^2(k), COBE-normalized, based on the approximations of
403Eisenstein & Hu (1997). k is in h^-1 Mpc. */
404/* TFmdm_set_cosm() and cobenorm() must have been called before this */
405/* Returns the density-weighted CDM+Baryon+Neutrino power spectrum */
406{
407 extern double tf_cbnu, growth_to_z0;
408 TFmdm_onek_hmpc(k); /* This sets the value in tf_cbnu */
409 return pow(2997.0*k, tilt+3.0)*SQR(deltaH*tf_cbnu*growth_to_z0);
410}
411
412double powerspcb(double k)
413/* Returns Delta^2(k), COBE-normalized, based on the approximations of
414Eisenstein & Hu (1997). k is in h^-1 Mpc. */
415/* TFmdm_set_cosm() and cobenorm() must have been called before this */
416/* Returns the density-weighted CDM+Baryon power spectrum */
417{
418 extern double growth_to_z0;
419 return pow(2997.0*k, tilt+3.0)*SQR(deltaH*TFmdm_onek_hmpc(k)*growth_to_z0);
420}
421
422double cobenorm()
423/* Return the Bunn & White (1997) fit for delta_H */
424/* Given lambda, omega_m, qtensors, and tilt */
425/* Open model with tensors is from Hu & White */
426{
427 double n;
428 n = tilt-1;
429 if (fabs(omega_m+lambda-1.0)<1e-5) { /* Flat universe */
430 if (qtensors)
431 return 1.94e-5*pow(omega_m, -0.785-0.05*log(omega_m))*
432 exp(n+1.97*n*n);
433 else return 1.94e-5*pow(omega_m, -0.785-0.05*log(omega_m))*
434 exp(-0.95*n-0.169*n*n);
435 /* Error -- was 0.0169 in earlier version of code */
436 } else if (fabs(lambda)<1e-5) { /* No lambda */
437 if (qtensors)
438 return 1.95e-5*pow(omega_m,-0.35-0.19*log(omega_m)-0.15*n)*
439 exp(+1.02*n+1.7*n*n);
440 /* Error -- was exp(-1.02n-1.7*n*n) in earlier version (v1-1.3)*/
441 else return 1.95e-5*pow(omega_m, -0.35-0.19*log(omega_m)-0.17*n)*
442 exp(-n-0.14*n*n);
443 /* Error -- was exp(-n-0.14*n*n) in earlier version (v1-1.2)*/
444 /* Error -- was exp(n+0.14*n*n) in earlier version (v1.3) */
445 } else return 1e-5*(2.422-1.166*exp(omega_m)+0.800*exp(lambda)
446 +3.780*omega_m-2.267*omega_m*exp(lambda)+0.487*SQR(omega_m)+
447 0.561*lambda+3.329*lambda*exp(omega_m)-8.568*omega_m*lambda+
448 1.080*SQR(lambda));
449}
450
451/* The following routine is from Numerical Recipes in C, 2nd edition,
452by Press, Teulkolsky, Vetterling, and Flannery, (Cambridge Univ. Press) */
453double erfcc(double x)
454{
455 double t,z,ans;
456 z=fabs(x);
457 t=1.0/(1.0+0.5*z);
458 ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
459 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
460 t*(-0.82215223+t*0.17087277)))))))));
461 return x >= 0.0 ? ans : 2.0-ans;
462}
463
464/* ======================= Integration Engine ======================= */
465
466double midpoint_int(double (*func)(double))
467/* Integrate the function func() on the range (0,1] */
468/* Use midpoint rule, expanding by factor of 3 each time */
469/* Linearly extrapolate to zero stepsize. Not fancy */
470{
471 int iter, bins, j;
472 double width, xeval, sum[11];
473
474 sum[0] = (*func)(0.5); /* First step is easy */
475 bins = 1; width = 1.0;
476 for (iter=1; iter<11; iter++) { /* Now start increasing the bins */
477 bins *= 3;
478 width /= 3.0;
479 sum[iter]=0.0;
480 for (xeval=width/2.0,j=1;j<=bins;xeval+=width,j++) {
481 if (j%3==2) continue; /* We've already done this one */
482 sum[iter]+=(*func)(xeval);
483 }
484 /* Add the new evaluations to the old ones. */
485 sum[iter]=sum[iter]*width+sum[iter-1]/3.0;
486 /* printf("%d iteration, sum = %f\n", iter, sum[iter]); */
487 /* Now decide if we're done */
488 if (iter<4) continue; /* We certainly want to do some work */
489 /* If error scales as width^2, then linear extrapolation in
490 quantity width^2 is sum[iter]+(sum[iter]-sum[iter-1])/8.0 */
491 if (fabs( (sum[iter]+(sum[iter]-sum[iter-1])/8.0)/
492 (sum[iter-1]+(sum[iter-1]-sum[iter-2])/8.0)-1.0)<1e-4)
493 return sum[iter]+(sum[iter]-sum[iter-1])/8.0;
494 /* if (fabs(sum[iter-1]/sum[iter]-1)<1e-4) return sum[iter]; */
495 }
496 crash("Too many steps in the integration.");
497}
498
499double integrate(double (*func)(double))
500/* Integrate the function func() on the range (0,1] */
501{
502 double midpoint_int(double (*func)(double));
503 return midpoint_int(func);
504}
505
506/* The above is a simple integrator. You could get faster execution times
507by using a more accurate integrator. The following is an example from
508Numerical Recipes, but all that is needed is for the call to integrate()
509to return the value of func() integrated on the range (0,1]. */
510
511#ifdef YES_NR
512#include "qromo.c"
513#include "midpnt.c"
514#include "polint.c"
515#include "nrutil.c"
516
517float integrate(float (*func)(float))
518/* Integrate the function func() on the range (0,1] */
519{
520 float qromo(float (*func)(float), float a, float b,
521 float (*choose)(float(*)(float),float,float,int));
522 float midpnt(float (*func)(float), float a, float b, int n);
523 return qromo(func,0,1,midpnt);
524}
525#endif
526
Note: See TracBrowser for help on using the repository browser.