From 269093d75e415a5ffcbae4f564f0e1363a15280a Mon Sep 17 00:00:00 2001 From: higginbottom Date: Thu, 27 Jul 2017 10:59:06 +0100 Subject: [PATCH 1/8] Fixed emittance_bb - there is now a check to see if we are about to integrate past ALPHA=100 --- source/bb.c | 11 +++++++++-- source/get_atomicdata.c | 4 +--- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/source/bb.c b/source/bb.c index 2b4071590..4e5730a94 100755 --- a/source/bb.c +++ b/source/bb.c @@ -625,11 +625,18 @@ emittance_bb (freqmin, freqmax, t) - if (alphamin > ALPHAMIN && alphamax < ALPHAMAX) + if (alphamin > ALPHAMIN && alphamax < ALPHAMAX) //We are within the tabulated range { return (q1 * t * t * t * t * integ_planck_d (alphamin, alphamax)); } - else + else if (alphamax > ALPHABIG) + { + if (alphamin > ALPHABIG) //The whole band is above the poitn where we can sensibly integrate the BB function + return(0); + else //only the upper part of the band is above ALPHABIG + return (q1 * t * t * t * t * qromb (planck_d, alphamin, ALPHABIG, 1e-7)); + } + else //We are outside the tabulated range and must integrate { return (q1 * t * t * t * t * qromb (planck_d, alphamin, alphamax, 1e-7)); } diff --git a/source/get_atomicdata.c b/source/get_atomicdata.c index 543218c8d..d66f571b3 100755 --- a/source/get_atomicdata.c +++ b/source/get_atomicdata.c @@ -2979,9 +2979,7 @@ index_phot_top () for (n = 0; n < ntop_phot + nxphot; n++) { - phot_top_ptr[n] = &phot_top[index[n + 1] - 1]; - printf ("NPHOT test n=%i f_0=%e z=%i state=%i\n",n,phot_top_ptr[n]->freq[0],phot_top_ptr[n]->z,phot_top_ptr[n]->istate); - + phot_top_ptr[n] = &phot_top[index[n + 1] - 1]; } /* Free the memory for the arrays */ From 99b6d9e8c4cb5c349e38874816c4d5e00d43b2e9 Mon Sep 17 00:00:00 2001 From: higginbottom Date: Sat, 29 Jul 2017 07:14:14 +0100 Subject: [PATCH 2/8] to test --- source/emission.c | 28 +++++++++++++++++++++++++++- source/py_wind_sub.c | 14 +++++++++++--- source/recomb.c | 12 +++++++++--- 3 files changed, 47 insertions(+), 7 deletions(-) diff --git a/source/emission.c b/source/emission.c index b3d47a0df..e62ca7428 100755 --- a/source/emission.c +++ b/source/emission.c @@ -199,13 +199,31 @@ total_emission (one, f1, f2) integrated */ { double t_e; - int nplasma; + int nplasma,n; PlasmaPtr xplasma; nplasma = one->nplasma; xplasma = &plasmamain[nplasma]; t_e = xplasma->t_e; // Change so calls to total emission are simpler + + + +//f2=1e16; +//f1=1e15; +//t_e=400; + +//for (n=200;n<600;n++) +//{ +// t_e=pow(10.,n/100.); +// printf ("FB_FULL temp %e total_fb %e %e\n",t_e,total_fb (one, t_e, f1, f2, FB_FULL, OUTER_SHELL),total_fb (one, t_e, f1, f2, FB_REDUCED, OUTER_SHELL)); + + //} + //exit(0); + + + + if (f2 < f1) { @@ -391,14 +409,20 @@ for (n=0;nt_e=5000; + + //!BUG SSMay04 //It doesn't seem to work unless this is zero? (SS May04) @@ -653,10 +659,10 @@ use that instead if possible -- 57h */ /* At this point, the variable nnn stores the number of points */ +// for (n=0;nt_e,fb_x[n],fb_y[n]); - - - +// exit(0); if (cdf_gen_from_array (&cdf_fb, fb_x, fb_y, nnn, f1, f2) != 0) { From 90a7f9afd4ae9fdb52da1b594aa1c90b2f4723ec Mon Sep 17 00:00:00 2001 From: higginbottom Date: Tue, 1 Aug 2017 10:00:57 +0100 Subject: [PATCH 3/8] more work --- source/cdf.c | 6 +++--- source/emission.c | 12 ++++++------ source/py_wind_sub.c | 6 ++++++ source/wind_updates2d.c | 3 ++- 4 files changed, 17 insertions(+), 10 deletions(-) diff --git a/source/cdf.c b/source/cdf.c index 44abbb164..f3cd1a9ba 100755 --- a/source/cdf.c +++ b/source/cdf.c @@ -371,9 +371,9 @@ cdf_gen_from_func (cdf, func, xmin, xmax, njumps, jump) } // 57ib - printf ("BLAH %e %e %e\n",cdf->x[0],cdf->y[0],cdf->d[0]); - printf ("BLAH %e %e %e\n",cdf->x[1],cdf->y[1],cdf->d[1]); - printf ("BLAH %e %e %e\n",cdf->x[2],cdf->y[2],cdf->d[2]); +// printf ("BLAH %e %e %e\n",cdf->x[0],cdf->y[0],cdf->d[0]); +// printf ("BLAH %e %e %e\n",cdf->x[1],cdf->y[1],cdf->d[1]); +// printf ("BLAH %e %e %e\n",cdf->x[2],cdf->y[2],cdf->d[2]); /* Check the pdf */ diff --git a/source/emission.c b/source/emission.c index e62ca7428..3056750e4 100755 --- a/source/emission.c +++ b/source/emission.c @@ -409,20 +409,20 @@ for (n=0;nt_e, xplasma->lum_tot_ioniz, xplasma->lum_lines_ioniz, xplasma->lum_ff_ioniz, xplasma->lum_rr_ioniz); + Log + ("t_e %8.2e lum_tot %8.2e lum_lines %8.2e lum_ff %8.2e lum_rr %8.2e \n", + xplasma->t_e, xplasma->lum_tot, xplasma->lum_lines, xplasma->lum_ff, xplasma->lum_rr); Log ("t_e %8.2e cool_tot %8.2e lum_lines %8.2e lum_ff %8.2e cool_rr %8.2e cool_comp %8.2e cool_adiab %8.2e cool_DR %8.2e \n", xplasma->t_e, xplasma->cool_tot_ioniz + xplasma->cool_comp_ioniz + xplasma->cool_adiabatic_ioniz + xplasma->cool_dr_ioniz, diff --git a/source/wind_updates2d.c b/source/wind_updates2d.c index 38a4a5cad..caf14d82a 100644 --- a/source/wind_updates2d.c +++ b/source/wind_updates2d.c @@ -725,7 +725,8 @@ WindPtr (w); /* JM130621- bugfix for windsave bug- needed so that we have the luminosities from ionization cycles in the windsavefile even if the spectral cycles are run */ - + if (nplasma==62) + printf("BLAH nwing=%i nplasma=%i lum_rr=%e\n",plasmamain[nplasma].nwind,nplasma,plasmamain[nplasma].lum_rr); plasmamain[nplasma].cool_tot_ioniz = plasmamain[nplasma].cool_tot; plasmamain[nplasma].lum_ff_ioniz = plasmamain[nplasma].lum_ff; plasmamain[nplasma].cool_rr_ioniz = plasmamain[nplasma].cool_rr; From 20f4e95e5a887107976bfac80ceaced10299a86d Mon Sep 17 00:00:00 2001 From: higginbottom Date: Tue, 1 Aug 2017 10:36:15 +0100 Subject: [PATCH 4/8] made changes to fix bug #284 --- source/ionization.c | 1 + source/wind_updates2d.c | 12 ++++++++---- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/source/ionization.c b/source/ionization.c index cd5f80f23..d1ec2cedb 100644 --- a/source/ionization.c +++ b/source/ionization.c @@ -430,6 +430,7 @@ one_shot (xplasma, mode) { xplasma->t_e = TMAX; } + zero_emit(xplasma->t_e); //Get the heating and cooling rates correctly for the new temperature } diff --git a/source/wind_updates2d.c b/source/wind_updates2d.c index caf14d82a..45e596e4c 100644 --- a/source/wind_updates2d.c +++ b/source/wind_updates2d.c @@ -689,6 +689,12 @@ WindPtr (w); } /* Check the balance between the absorbed and the emitted flux */ + + //NSH 0717 - first we need to ensure the cooling and luminosities reflect the current temperature + + cool_sum = wind_cooling (0.0, VERY_BIG); /*We call wind_cooling here to obtain an up to date set of cooling rates */ + lum_sum = wind_luminosity (0.0, VERY_BIG); /*and we also call wind_luminosity to get the luminosities */ + xsum = psum = ausum = lsum = fsum = csum = icsum = apsum = aausum = abstot = 0; //1108 NSH zero the new csum counter for compton heating @@ -725,8 +731,7 @@ WindPtr (w); /* JM130621- bugfix for windsave bug- needed so that we have the luminosities from ionization cycles in the windsavefile even if the spectral cycles are run */ - if (nplasma==62) - printf("BLAH nwing=%i nplasma=%i lum_rr=%e\n",plasmamain[nplasma].nwind,nplasma,plasmamain[nplasma].lum_rr); + plasmamain[nplasma].cool_tot_ioniz = plasmamain[nplasma].cool_tot; plasmamain[nplasma].lum_ff_ioniz = plasmamain[nplasma].lum_ff; plasmamain[nplasma].cool_rr_ioniz = plasmamain[nplasma].cool_rr; @@ -776,8 +781,7 @@ WindPtr (w); - cool_sum = wind_cooling (0.0, VERY_BIG); /*We call wind_cooling here to obtain an up to date set of cooling rates */ - lum_sum = wind_luminosity (0.0, VERY_BIG); /*and we also call wind_luminosity to get the luminosities */ + if (modes.zeus_connect == 1 && geo.hydro_domain_number > -1) //If we are running in zeus connect mode, we output heating and cooling rates. From 05f882f6f4bf1eeb8909172a1714dc42a613e4c9 Mon Sep 17 00:00:00 2001 From: higginbottom Date: Wed, 2 Aug 2017 11:51:55 +0100 Subject: [PATCH 5/8] testing --- source/cdf.c | 15 ++++++++++++--- source/continuum.c | 46 ++++++++++++++++++++++++++++++++++++++++---- source/py_wind_sub.c | 10 +++++----- 3 files changed, 59 insertions(+), 12 deletions(-) diff --git a/source/cdf.c b/source/cdf.c index f3cd1a9ba..d77bf5a14 100755 --- a/source/cdf.c +++ b/source/cdf.c @@ -569,7 +569,7 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax) double xmin, xmax; { int allzero; - int m, n, nn; + int m, n; double sum, q; int echeck,cdf_check (), recalc_pdf_from_cdf (); @@ -624,7 +624,7 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax) //Now at the end - nn = n_xy; + n = n_xy-1; while (y[n] == 0.0) { xmax = x[n]; @@ -718,7 +718,10 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax) sum = pdf_z[pdf_n - 1]; for (n = 1; n < pdf_n; n++) + { + printf ("%e %e\n",pdf_x[n],pdf_z[n]); pdf_z[n] /= sum; + } @@ -781,7 +784,13 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax) } // 57ib if ((echeck = cdf_check (cdf)) != 0) { - Error ("cdf_gen_from_array: error %d on cdf_check\n", echeck); + Error ("cdf_gen_from_array: error %d on cdf_check\n", echeck); + for (n=0;n lambdamin && comp[spectype].xmod.w[n] <= lambdamax) + { + w_local[nwave]=comp[spectype].xmod.w[n]; + f_local[nwave]=comp[spectype].xmod.f[n]; + nwave++; + } + } + + if (comp[spectype].xmod.w[0] < lambdamax && lambdamax < comp[spectype].xmod.w[comp[spectype].nwaves-1]) + { + w_local[nwave]=lambdamax; + linterp(lambdamax, comp[spectype].xmod.w, comp[spectype].xmod.f, comp[spectype].nwaves, &y, 0); + f_local[nwave]=y; + nwave++; + } +// printf ("BLAH min %e max %e\n",lambdamin,lambdamax); +// printf ("BLAH wmin %e wmax %e\n",comp[spectype].xmod.w[0],comp[spectype].xmod.w[comp[spectype].nwaves-1]); + +// for (n=0;n Date: Thu, 3 Aug 2017 11:37:10 +0100 Subject: [PATCH 6/8] Fixed comments in cdf.c --- source/cdf.c | 100 +++++++++++++++++++++++++++------------------------ 1 file changed, 53 insertions(+), 47 deletions(-) diff --git a/source/cdf.c b/source/cdf.c index 54a4d5478..e6fb5c0f5 100755 --- a/source/cdf.c +++ b/source/cdf.c @@ -8,14 +8,14 @@ The main routines are - pdf_gen_from_func (&pdf,&func,xmin,xmax,njumps,jump) + cdf_gen_from_func (&cdf,&func,xmin,xmax,njumps,jump) Generate a cumulative distrubution function (cdf) from a function - pdf_gen_from_array(&pdf,x,y,n_xy,xmin,xmax,njumps,jump) + cdf_gen_from_array(&cdf,x,y,n_xy,xmin,xmax,njumps,jump) Generate a cumulative distribution function (cdf) from an array of values - pdf_get_rand(&pdf) + cdf_get_rand(&cdf) Generate a single sample from a cdf defined by either of the first two routines @@ -27,10 +27,10 @@ pdf. Note that pdf_gen-from_func or pdf_gen_from_array should have been called previously! - pdf_limit(&pdf,xmin,xmax) sets limit1 and limit2 so that one can generate - distributions from a limited range within a previously created pdf. + cdf_limit(&pdf,xmin,xmax) sets limit1 and limit2 so that one can generate + distributions from a limited range within a previously created cdf. - pdf_get_rand_limit(&pdf) will samples the pdf but limits the range of + cdf_get_rand_limit(&pdf) will samples the cdf but limits the range of the return to lie between xmin and xmax from pdf_limit. @@ -41,18 +41,18 @@ number of points. There is also a simple diagnostic routine - pdf_check(&pdf) + cdf_check(&cdf) - which returns a nonzero number if certain checks of the pdf fail + which returns a nonzero number if certain checks of the cdf fail and a routine - pdf_to_file(&pdf,filename) - to write a pdf structure to a file + cdf_to_file(&cdf,filename) + to write a cdf structure to a file -Arguments for pdf_gen_from_func - PdfPtr pdf; A ptr to a pdf structure (see below) +Arguments for cdf_gen_from_func + CdfPtr cdf; A ptr to a cdf structure (see below) double (*func)() The probablility density function to be integrated to create the pdf double xmin,xmax; The range over which the function is to be @@ -63,32 +63,30 @@ Arguments for pdf_gen_from_func and xmax int njumps; The size of jumps[] -Arguments for pdf_gen_from_array are same as above except +Arguments for cdf_gen_from_array are same as above except - double x[],y[] The one dimensional arrays containg the places x where the cumulative - distribution function density has been evaluated. The array need not - be uniformly spaced. y is the value of the cdf at the point x - double d[] The rate of change of the probability density at x, which is calculated - from the CDF. -- Added 57i + double x[],y[] The one dimensional arrays containg the places x where the probability + density function density has been evaluated. The array need not + be uniformly spaced. y is the value of the pdf at the point x int n_xy The size of x[] and y[] Arguments for pdf_get_rand - PdfPtr pdf pdf (See below) is a structure which contains the cumulative + CdfPtr cdf cdf (See below) is a structure which contains the cumulative distribution function. -Arguments for pdf_limit - double xmin,xmax Here xmin and xmax will cause pdf[].limit1 and pdf[].limit2 to be set +Arguments for cdf_limit + double xmin,xmax Here xmin and xmax will cause cdf[].limit1 and cdf[].limit2 to be set in such a way that when you call pdf_get_rand_limit the distribution will be sampled only between xmin and xmax. Watch out though because if xmin and xmax lie outside of the original value of xmin and xmax then the distribution will not extend this far. Arguments for pdf_to_file - PdfPtr pdf + PdfPtr cdf char filename[] If you can't figure this out, you are in deep trouble! Returns: - All of the routines return 0 on successful completion, except pdf_get_rand which returns - a random value between xmin and xmax. Multiple calls to pdf_get_rand should regenerate + All of the routines return 0 on successful completion, except cdf_get_rand which returns + a random value between xmin and xmax. Multiple calls to cdf_get_rand should regenerate the probability density function. If there is an error the routines usually send an error message to the screen and logfile @@ -101,38 +99,36 @@ Arguments for pdf_to_file and then allow one to sample that distribution function. The procedure is either to call pdf_gen_from_func if you have a function whose probability distribution you want (as might be the case if you were generating bremsstrahlung photons and you knew - the formula) or pdf_gen_from_array (as might be the case if you were trying to + the formula) or cdf_gen_from_array (as might be the case if you were trying to generate Monte Carlo spectra from precalculated Kurucz models). The result of calling either of these routines will be to populate a structure of - the form Pdf - pdf->x[] will contain NPDF values of x between xmin and xmax. The values will be + the form Cdf + cdf->x[] will contain NPDF values of x between xmin and xmax. The values will be spaced so that it is almost but not quite true that if you generate a uniform number n between 0 and NPDF you could just read off the value pdf[n].x and that by doing this multiple times you could regenerate the distribution function - pdf->y[] will contain values between 0 and 1 and is the integral of the cumulative + cdf->y[] will contain values between 0 and 1 and is the integral of the cumulative distribution function from xmin to pdf[].x - Once a pdf has been generated, one samples the pdf by calls to pdf_get_rand which - creates a random number between 0 and 1, finds the elements in pdf which surround + Once a cdf has been generated, one samples the cdf by calls to cdf_get_rand which + creates a random number between 0 and 1, finds the elements in cdf which surround the random number and interpolates to return a value x. It is possible to force specific values of x to appear in the pdf. This is desirable if there are edges where the probability density changes rapidly, as for example at the Lyman edge in a stellar spectrum. By using jumps properly you can assure - that edges will be sharp in the regenerated distribution (since otherwse pdf_get_rand + that edges will be sharp in the regenerated distribution (since otherwse cdf_get_rand will linearly interpolate between two points on the opposite sides of the edge). Notes: Must be compiled with log.c since it makes use of the "python" logging routines for errors. - The fact that these routines start with pdf... is historical and represents a poor - choice of nomenclature. We reallly mean comulative distribution functions - It would probably be desirable to modify this so that the size of the parallel arrays being generated could be determined from inside the program. In that case - NPDF would have to become the maximum size. + NPDF would have to become the maximum size. - NSH stis was tried in Jul 2017 - turns + out to be very hard! Error??: For Kurucz models, The probability distribution beyond the He edge @@ -142,7 +138,9 @@ Arguments for pdf_to_file everything down. This is an another argument that the whole routine pdf_gen_from_array neeeds to be rewritten to avoid using NPDF altogether. One way to do this would be simply to use the values of x and y in creating an inital pdf. There - is however a question of resampling. + is however a question of resampling. - NSH Jul 2017 - this has been addressed + by removing the sampling aspect. Turns out with the faster machines we dont + need to be so careful about making it easy to find the right point in a CDF. History: @@ -180,6 +178,9 @@ Arguments for pdf_to_file the number of points that are used in the array pdf_array in situations where the binning is too course. In the process, eliminated pdf_init. It was only called by one routine. + 17Jul nsh Modified to make everything refer to CDFs as they should be, also changed + cdf_gen_from_array to work on a supplied array - no jumps. + **************************************************************/ @@ -215,7 +216,7 @@ double *pdf_array; Description: -This routine stores values of the function in pdf_array +This routine stores values of the function in cdf_array 02feb ksl Modified to handle a problem in which the pdf was not monotonic because of the @@ -391,7 +392,7 @@ cdf_gen_from_func (cdf, func, xmin, xmax, njumps, jump) gen_array_from_func -This is a routine which is called by pdf_gen_from_func which simply calculates the cumulative +This is a routine which is called by cdf_gen_from_func which simply calculates the cumulative distribution of the function in equally spaced steps between xmin and xmax. The CDF is properly normalized. @@ -515,7 +516,7 @@ gen_array_from_func (func, xmin, xmax, pdfsteps) /* -pdf_gen_from_array +cdf_gen_from_array The one dimensional arrays containg the places x where the probability density y has been evaluated. The array need not be uniformly spaced. @@ -525,7 +526,7 @@ one can linearly interpolate between points Notes: - 06sep -- There are differences between this and pdf_gen_from_func + 06sep -- There are differences between this and cdf_gen_from_func that look historical -- ksl @@ -553,7 +554,7 @@ one can linearly interpolate between points so that problems with this would be easier to update in future. 1405 JM -- Increased PDF array for use with disk14 models - 1707 NSH -- Removed code for jumps - we now just supply asuitable unscaled pdf + 1707 NSH -- Removed code for jumps - we now just supply a suitable unscaled pdf */ #define PDF_ARRAY 28000 @@ -801,7 +802,7 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax) /* -pdf_get_rand +cdf_get_rand Our searching routine assumes we can predict roughly where in the @@ -812,7 +813,8 @@ structure the right values will be 06sep ksl 57i -- Modified to account for the fact that the probability density is not uniform between intervals. - 17jun ksl Modified for variable sizes of distribution function + 17jun ksl Modified for variable sizes of distribution function + 17jul nsh Modified to use gsl routine to find the right point in the CDF. */ double @@ -876,9 +878,8 @@ find that minimum and maximum value. Having done this one can call pdf_get_rand_ Error Some of this is not quite right. It's intended to let us generate limits in the CDF but we need x too - 06sep ksl 57i -- Added the xlimits to the pdf + 06sep ksl 57i -- Added the xlimits to the cdf 17jun ksl Modified for variable sizes of distribution function - 7 Jul nsh Modified to make everything refer to CDFs as they should be */ int @@ -954,7 +955,6 @@ History that the probability density is not uniform between intervals. 17jul ksl Modified for variable lengths of pdfs - 17jul nsh Changed terminology to CDF */ double @@ -1158,7 +1158,13 @@ cdf_check (cdf) allow one to calculate a first order correction to a uniform distibution between the points where the gradient has been calculated. - + + XXX NSH - the way that the ends are dealt with is bot + great - we should really try to come up with an + extrapolation rather than just fill in the same gradients + for the second and penultimate cells into the + forst and last cells. + Notes: d is calculated by differencing the cdf. Thus From 37158e323f4dd418e73d8aae29bd0bd81242472a Mon Sep 17 00:00:00 2001 From: higginbottom Date: Fri, 4 Aug 2017 17:10:38 +0100 Subject: [PATCH 7/8] fixed continuum CDF --- source/cdf.c | 224 ++++++++++++++++--------------------------- source/continuum.c | 28 +++++- source/cooling.c | 18 ++-- source/emission.c | 8 +- source/py_wind_sub.c | 9 +- source/recomb.c | 12 ++- 6 files changed, 137 insertions(+), 162 deletions(-) diff --git a/source/cdf.c b/source/cdf.c index e6fb5c0f5..63e45b9c3 100755 --- a/source/cdf.c +++ b/source/cdf.c @@ -570,9 +570,14 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax) double xmin, xmax; { int allzero; + int nmin,nmax,cdf_n; int m, n; double sum, q; - int echeck,cdf_check (), recalc_pdf_from_cdf (); +// int echeck,cdf_check (), recalc_pdf_from_cdf (); + int echeck; + + + if (n_xy > NCDF) @@ -590,6 +595,16 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax) } allzero = 0; + + if (x[0]!=xmin || x[n_xy-1]!=xmax) + { + Error("cdf_gen_from_array: input array does not run exactly from xmin to xmax %e != %e or %e != %e\n",x[0],xmin,x[n_xy-1],xmax); + for (n=0;n x[n_xy - 1] || allzero == 0) - { // These are special (probably nonsensical) cases - pdf_x[0] = xmin; - pdf_z[0] = 0.; - pdf_x[1] = xmax; - pdf_z[0] = 1.; - sum = 1.0; - pdf_n = 2; - Error - ("pdf_gen_from_array: all y's were zero or xmin xmax out of range of array x-- returning uniform distribution %d\n", - allzero); + n = n_xy-1; + nmax=n_xy; + while (y[n] == 0.0) + { + nmax = n; + n--; + } - } - else - { - m = 0; - while (x[m] <= xmin) - m++; // Find the bottom boundary - pdf_x[0] = xmin; - if (m == 0) //This is hit if the bottom boundary is outside the supplied data - { - pdf_y[0] = y[0]; //Assume prob. density is constant outside array lims. - } - else //We interpolate between supplied data points to get the value of the PDF at exactly xmin - { - q = (xmin - x[m - 1]) / (x[m] - x[m - 1]); - pdf_y[0] = q * y[m] + (1. - q) * y[m - 1]; - } - pdf_n = 1; - // Completed first element; now do those that are completely in the grid - while (x[m] < xmax && m < n_xy) - { - pdf_x[pdf_n] = x[m]; - pdf_y[pdf_n] = y[m]; - m++; - pdf_n++; - if (pdf_n > PDF_ARRAY) - { - Error - ("cdf_gen_from_array: pdf_n (%d) exceeded maximum array size PDF_ARRAY (%d) \n", - pdf_n, PDF_ARRAY); - Error - ("cdf_gen_from_array: n_xy %d xmin %f xmax %f\n", - n_xy, xmin, xmax); - Error - ("cdf_gen_from_array: Consider increasing PDF_ARRAY to a value > n_xy + njumps, and recompiling\n"); - exit (0); - } - } - // Now worry about the last element - pdf_x[pdf_n] = xmax; - m--; /* Reduce m so that x[m] is les than xmax ksl 170801 */ - if (m < n_xy - 1) //We are not at the limit of the supplied data, so we interpolate - { - q = (xmax - x[m]) / (x[m + 1] - x[m]); - pdf_y[pdf_n] = q * y[m + 1] + (1. - q) * y[m]; - } - else //xmax is outside the supplied data +nmax--; + + if (nmax==nmin) { - pdf_y[pdf_n] = y[m - 1]; //Again assume constant prob. density outside lims + Error ("cdf_gen_from_array - only one point in supplied PDF\n"); + exit(0); } - pdf_n++; - - + if (xmax < x[0] || xmin > x[n_xy - 1] || allzero == 0 ) + { // These are special (probably nonsensical) cases + cdf->x[0] = xmin; + cdf->y[0] = 0.; + cdf->x[1] = xmax; + cdf->y[1] = 1.; + cdf->norm = 1.0; + cdf->ncdf = 1; + Error + ("cdf_gen_from_array: all y's were zero or xmin xmax out of range of array x-- returning uniform distribution %d\n", + allzero); + } + else + { /* So at this point, have unscaled probability density in pdf_x, pdf_y for the points * specified by the input array but we want the cumulative distribution * We also have assured that there is one value of pdf_x that corresponds to all of the jumps */ + + // The following lines perform an integration via the trapezoid rule - each point contains + // the integral up to that poont, so it starts at 0 and ends at the total + + cdf_n=(nmax-nmin); - pdf_z[0] = 0.; - for (n = 1; n < pdf_n; n++) + cdf->x[0] = x[nmin]; + cdf->y[0] = 0.0; + for (n = 1; n < cdf_n+1; n++) { - pdf_z[n] = - pdf_z[n - 1] + 0.5 * (pdf_y[n - 1] + pdf_y[n]) * (pdf_x[n] - - pdf_x[n - 1]); + cdf->x[n]=x[nmin+n]; + cdf->y[n] = + cdf->y[n-1] + 0.5 * (y[nmin + n - 1] + y[nmin + n]) * (x[nmin + n] - + x[nmin + n - 1]); } - sum = pdf_z[pdf_n - 1]; + + sum = cdf->y[cdf_n - 1]; //the total integrated pdf - for (n = 1; n < pdf_n; n++) + for (n = 1; n < cdf_n; n++) { - printf ("%e %e\n",pdf_x[n],pdf_z[n]); - pdf_z[n] /= sum; + cdf->y[n] /= sum; //this is now a cdf - we go from 0 to 1. } - - -/* So pdf_z contains a properly normalized cdf on the points specified by - the input array, or more explicitly, at the points specied in the array - pdf_x -*/ - - - /* Add a check that the pdf_z is monotonic. This check should not really be necessary - * since by construction this should be the case*/ - - echeck = 0; - for (n = 1; n < pdf_n; n++) - { - if (pdf_z[n] < pdf_z[n - 1]) - { - Error ("cdf_gen_from_array: pdf_z is not monotonic at %d\n", n); - echeck = 1; - } - } - - if (echeck) - { - for (n = 0; n < pdf_n; n++) - { - Log ("cdf_gen_from_array: %5d %11.6e %11.6e %11.6e\n", n, - pdf_x[n], pdf_y[n], pdf_z[n]); - } - echeck = 0; - } - - } - - - - //We now simply put the CDF into the CDF array - - cdf->x[0] = xmin; - cdf->y[0] = 0; - for (n=1;nx[n] = pdf_x[n]; - cdf->y[n] = pdf_z[n]; - } - cdf->ncdf = pdf_n; - cdf->x[n] = xmax; - cdf->y[n] = 1.0; - cdf->norm = sum; - + cdf->ncdf = cdf_n; + cdf->x[cdf->ncdf] = x[nmax]; + cdf->y[cdf->ncdf] = 1.0; + cdf->norm = sum; + +} @@ -789,7 +725,8 @@ cdf_gen_from_array (cdf, x, y, n_xy, xmin, xmax) Error ("cdf_gen_from_array: error %d on cdf_check\n", echeck); for (n=0;nncdf+1;n++) + printf ("cdf_n=%i %e %e\n",n,cdf->x[n],cdf->y[n]); exit(0); @@ -1067,6 +1004,7 @@ cdf_check (cdf) ("cdf_check: cumulative distribution function should start at 0 not %e\n", y); hcheck = 1; + printf ("%e %e %e %e\n",cdf->x[0],cdf->y[0],cdf->x[1],cdf->y[1]); } if (cdf->y[cdf->ncdf] != 1.0) { diff --git a/source/continuum.c b/source/continuum.c index 8903c8bb1..2c3380178 100755 --- a/source/continuum.c +++ b/source/continuum.c @@ -89,13 +89,19 @@ one_continuum (spectype, t, g, freqmin, freqmax) for (n=0;n lambdamin && comp[spectype].xmod.w[n] <= lambdamax) { +// printf (" inserted %e %e\n",lambdamin,lambdamax); w_local[nwave]=comp[spectype].xmod.w[n]; f_local[nwave]=comp[spectype].xmod.f[n]; nwave++; } +// else +// { +// printf ("\n "); +// } + } if (comp[spectype].xmod.w[0] < lambdamax && lambdamax < comp[spectype].xmod.w[comp[spectype].nwaves-1]) @@ -105,6 +111,26 @@ one_continuum (spectype, t, g, freqmin, freqmax) f_local[nwave]=y; nwave++; } + + + //There are two pathological cases to deal with, when we only have one non zero point, we need to make an extra point just up/dpown from the penultimate/zecond point so we can make a sensible CDF. + + + if (f_local[nwave-2]==0.0) //We have a zero just inside the end + { + nwave++; + w_local[nwave-1]=w_local[nwave-2]; + f_local[nwave-1]=f_local[nwave-2]; + w_local[nwave-2]=w_local[nwave-3]/(1.-DELTA_V/(2.*C)); + linterp(w_local[nwave-2], comp[spectype].xmod.w, comp[spectype].xmod.f, comp[spectype].nwaves, &y, 0); + f_local[nwave-2]=y; + } + + + + + + // printf ("BLAH min %e max %e\n",lambdamin,lambdamax); // printf ("BLAH wmin %e wmax %e\n",comp[spectype].xmod.w[0],comp[spectype].xmod.w[comp[spectype].nwaves-1]); diff --git a/source/cooling.c b/source/cooling.c index 8c6f62d9c..dd4a14a35 100644 --- a/source/cooling.c +++ b/source/cooling.c @@ -153,8 +153,10 @@ xtotal_emission (one, f1, f2) { double t_e; int nplasma; + double cooling; PlasmaPtr xplasma; + cooling=0.0; nplasma = one->nplasma; xplasma = &plasmamain[nplasma]; @@ -172,36 +174,36 @@ xtotal_emission (one, f1, f2) //The first term here is the fb cooling due to macro ions and the second gives //the fb cooling due to simple ions. //total_fb has been modified to exclude recombinations treated using macro atoms. - xplasma->cool_tot = xplasma->cool_rr; + cooling = xplasma->cool_rr; //Note: This the fb_matom call makes no use of f1 or f2. They are passed for //now in case they should be used in the future. But they could //also be removed. // (SS) xplasma->lum_lines = total_bb_cooling (xplasma, t_e); - xplasma->cool_tot += xplasma->lum_lines; + cooling += xplasma->lum_lines; /* total_bb_cooling gives the total cooling rate due to bb transisions whether they are macro atoms or simple ions. */ xplasma->lum_ff = total_free (one, t_e, f1, f2); - xplasma->cool_tot += xplasma->lum_ff; + cooling += xplasma->lum_ff; } else //default (non-macro atoms) (SS) { /*The line cooling is equal to the line emission */ - xplasma->cool_tot = xplasma->lum_lines = total_line_emission (one, f1, f2); + cooling = xplasma->lum_lines = total_line_emission (one, f1, f2); /* The free free cooling is equal to the free free emission */ - xplasma->cool_tot += xplasma->lum_ff = total_free (one, t_e, f1, f2); + cooling += xplasma->lum_ff = total_free (one, t_e, f1, f2); /*The free bound cooling is equal to the recomb rate x the electron energy - the boinding energy - this is computed with the FB_REDUCED switch */ - xplasma->cool_tot += xplasma->cool_rr = total_fb (one, t_e, f1, f2, FB_REDUCED, OUTER_SHELL); //outer shell recombinations + cooling += xplasma->cool_rr = total_fb (one, t_e, f1, f2, FB_REDUCED, OUTER_SHELL); //outer shell recombinations } } - return (xplasma->cool_tot); + return (cooling); } @@ -334,7 +336,7 @@ wind_cooling (f1, f2) if (wmain[n].vol > 0.0) { nplasma = wmain[n].nplasma; - cool += x = xtotal_emission (&wmain[n], f1, f2); + cool += x = cooling (&plasmamain[nplasma],plasmamain[nplasma].t_e); lum_lines += plasmamain[nplasma].lum_lines; cool_rr += plasmamain[nplasma].cool_rr; lum_ff += plasmamain[nplasma].lum_ff; diff --git a/source/emission.c b/source/emission.c index 3056750e4..fab161ad6 100755 --- a/source/emission.c +++ b/source/emission.c @@ -821,13 +821,15 @@ one_ff (one, f1, f2) if (xplasma->t_e != one_ff_te || f1 != one_ff_f1 || f2 != one_ff_f2) { /* Generate a new pdf */ - - dfreq = (f2 - f1) / ARRAY_PDF-1; - for (n = 0; n < ARRAY_PDF; n++) + dfreq = (f2 - f1) / (ARRAY_PDF-1); + for (n = 0; n < ARRAY_PDF-1; n++) { ff_x[n] = f1 + dfreq * n; ff_y[n] = ff (one, xplasma->t_e, ff_x[n]); } + + ff_x[ARRAY_PDF-1] = f2; + ff_y[ARRAY_PDF-1] = ff (one, xplasma->t_e, ff_x[ARRAY_PDF-1]); diff --git a/source/py_wind_sub.c b/source/py_wind_sub.c index 7fa0eee90..e69dfb2cb 100755 --- a/source/py_wind_sub.c +++ b/source/py_wind_sub.c @@ -3304,10 +3304,11 @@ ionH1\tionH2\tionHe1\tionHe2\tionHe3\tionC3\tionC4\tionC5\tionN5\tionO6\tionSi4\ w[n].thetacen / RADIAN, vtot, w[n].v[0], w[n].v[1], w[n].v[2], w[n].dvds_ave, w[n].vol, plasmamain[np].rho, plasmamain[np].ne, plasmamain[np].t_e, plasmamain[np].t_r, plasmamain[np].ntot, plasmamain[np].w, plasmamain[np].ave_freq, plasmamain[np].ip, plasmamain[np].xi, plasmamain[np].converge_whole, plasmamain[np].converge_t_r, plasmamain[np].converge_t_e, plasmamain[np].converge_hc, - plasmamain[np].cool_tot , - plasmamain[np].lum_tot, plasmamain[np].lum_rr,plasmamain[np].cool_rr, plasmamain[np].lum_ff, plasmamain[np].lum_lines, plasmamain[np].cool_adiabatic, - plasmamain[np].cool_comp, plasmamain[np].cool_dr, plasmamain[np].heat_tot, plasmamain[np].heat_photo, plasmamain[np].heat_auger, - plasmamain[np].heat_lines, plasmamain[np].heat_ff, plasmamain[np].heat_comp, plasmamain[np].heat_ind_comp, h1den, h2den, he1den, + plasmamain[np].cool_tot_ioniz+ plasmamain[np].cool_comp_ioniz + plasmamain[np].cool_adiabatic_ioniz + plasmamain[np].cool_dr_ioniz, + plasmamain[np].lum_tot_ioniz, plasmamain[np].lum_rr_ioniz,plasmamain[np].cool_rr_ioniz, plasmamain[np].lum_ff_ioniz, plasmamain[np].lum_lines_ioniz, + plasmamain[np].cool_adiabatic_ioniz, plasmamain[np].cool_comp_ioniz, plasmamain[np].cool_dr_ioniz, plasmamain[np].heat_tot, plasmamain[np].heat_photo, + plasmamain[np].heat_auger, plasmamain[np].heat_lines, plasmamain[np].heat_ff, plasmamain[np].heat_comp, plasmamain[np].heat_ind_comp, + h1den, h2den, he1den, he2den, he3den, c3den, c4den, c5den, n5den, o6den, si4den); } else diff --git a/source/recomb.c b/source/recomb.c index bb17b9958..332c0dc67 100644 --- a/source/recomb.c +++ b/source/recomb.c @@ -620,13 +620,13 @@ use that instead if possible -- 57h */ nnn=0; //Zero the index for elements in the flux array nn=0; //Zero the index for elements in the jump array n=0; //Zero the counting element for equally spaced frequencies - dfreq = (f2 - f1) / ARRAY_PDF; //This is the frequency spacing for the equally spaced elements + dfreq = (f2 - f1) / (ARRAY_PDF-1); //This is the frequency spacing for the equally spaced elements while (n < (ARRAY_PDF) && nnn < NCDF) //We keep going until n=ARRAY_PDF-1, which will give the maximum required frequency { freq=f1 + dfreq * n; //The frequency of the array element we would make in the normal run of things if (freq > fb_jumps[nn] && nn < fb_njumps) //The element we were going to make has a frequency abouve the jump { - fb_x[nnn]=fb_jumps[nn]*(1.-DELTA_V/(2.*C)); //We make one frequency point 1km/s below the jump + fb_x[nnn]=fb_jumps[nn]*(1.-DELTA_V/(2.*C)); //We make one frequency point DELTA_V cm/s below the jump fb_y[nnn]=fb (xplasma, xplasma->t_e, fb_x[nnn], nions, FB_FULL); //And the flux for that point nnn=nnn+1; //increase the index of the created array fb_x[nnn]=fb_jumps[nn]*(1.+DELTA_V/(2*C)); //And one frequency point just above the jump @@ -636,7 +636,7 @@ use that instead if possible -- 57h */ } else //We haven't hit a jump { - if (freq > fb_x[nnn-1]) //Deal with the unusual case where the upper point in our 'jump' pair is above the next ragular point + if (freq > fb_x[nnn-1]) //Deal with the unusual case where the upper point in our 'jump' pair is above the next regular point { fb_x[nnn] = freq; //Set the next array element frequency fb_y[nnn] = fb (xplasma, xplasma->t_e, fb_x[nnn], nions, FB_FULL); //And the flux @@ -649,6 +649,12 @@ use that instead if possible -- 57h */ } } } + + //Ensure the last point lines up exatly with f2 + + fb_x[nnn-1]=f2; + fb_y[nnn-1]=fb (xplasma, xplasma->t_e, f2, nions, FB_FULL); + if (nnn > NCDF) { From 2d05dbc0d456bb821234a5e2f5ea77791f5c7932 Mon Sep 17 00:00:00 2001 From: higginbottom Date: Mon, 7 Aug 2017 11:57:51 +0100 Subject: [PATCH 8/8] Some more small changes --- source/bb.c | 5 +++-- source/cdf.c | 11 +++++++---- source/cooling.c | 28 ++++++---------------------- 3 files changed, 16 insertions(+), 28 deletions(-) diff --git a/source/bb.c b/source/bb.c index 4e5730a94..851b7f5c4 100755 --- a/source/bb.c +++ b/source/bb.c @@ -163,6 +163,7 @@ check_fmax(fmin,fmax,temp) This is a little helper routine written by NSH in Au 12nov ksl Removed some of the old Notes assocaited with this routine. See version earlier than 74 for these old notes. 17jul nsh - changed references to PDFs to CDFs + 17jul nsh - also changed how the BB cdf is generated - low alpha jumps were not really doing anything! **************************************************************/ #define ALPHAMIN 0.4 // Region below which we will use a low frequency approximation @@ -607,7 +608,7 @@ planck_d (alpha) // Calculate the emittance of a bb between freqmin and freqmax // Should integrate to sigma -//NSH - 17Jul - made change to if/else statement so that qromb used whenever bounds are not completely withing the tabulated bands. +//NSH - 17Jul - made change to if/else statement so that qromb used whenever bounds are not completely within the tabulated bands. double emittance_bb (freqmin, freqmax, t) double freqmin, freqmax, t; @@ -631,7 +632,7 @@ emittance_bb (freqmin, freqmax, t) } else if (alphamax > ALPHABIG) { - if (alphamin > ALPHABIG) //The whole band is above the poitn where we can sensibly integrate the BB function + if (alphamin > ALPHABIG) //The whole band is above the point where we can sensibly integrate the BB function return(0); else //only the upper part of the band is above ALPHABIG return (q1 * t * t * t * t * qromb (planck_d, alphamin, ALPHABIG, 1e-7)); diff --git a/source/cdf.c b/source/cdf.c index 63e45b9c3..9dd947e92 100755 --- a/source/cdf.c +++ b/source/cdf.c @@ -554,7 +554,10 @@ one can linearly interpolate between points so that problems with this would be easier to update in future. 1405 JM -- Increased PDF array for use with disk14 models - 1707 NSH -- Removed code for jumps - we now just supply a suitable unscaled pdf + 1707 NSH -- Removed code for jumps - we now just supply a suitable unscaled pdf, + The code expects data between xmin and xmax - it scan cope with + leading ro training zeros, but does require at least two non zero + points, otherwise it cannt make a CDF. */ #define PDF_ARRAY 28000 @@ -712,9 +715,6 @@ nmax--; } - - - /* Calculate the gradients */ if (calc_cdf_gradient (cdf)) { @@ -1158,6 +1158,9 @@ calc_cdf_gradient (cdf) /* Fill in the ends */ cdf->d[0] = cdf->d[1]; + + //XXX NSH - this routine does not treat the ends of the array very well - it just asumes the same gradient - we could do better... + //NSH - a couple of other ways of doing this - linear interpolation from the previous points, or just the gradient of the current interval - none works great. // cdf->d[0]=cdf->d[1]-((cdf->d[2]-cdf->d[1])/(cdf->x[2]-cdf->x[1]))*(cdf->x[1]-cdf->x[0]); diff --git a/source/cooling.c b/source/cooling.c index dd4a14a35..0cfe69782 100644 --- a/source/cooling.c +++ b/source/cooling.c @@ -67,26 +67,8 @@ cooling (xxxplasma, t) /*81c - nsh - we now treat DR cooling as a recombinational process - still unsure as to how to treat emission, so at the moment it remains here */ - xxxplasma->cool_dr = total_fb (&wmain[xxxplasma->nwind], t, 0, VERY_BIG, FB_REDUCED, INNER_SHELL); - -// f1=2.641e15; -// f2=1.523e16; - -// df=(f2-f1)/200; -// printf ("FB_TEST total %e\n",total_fb (&wmain[xxxplasma->nwind], t, f1, f2, FB_FULL, OUTER_SHELL)); -// for (n=1;n<200;n++) -// { -// lower= total_fb (&wmain[xxxplasma->nwind], t, f1, f1+(n*df), FB_FULL, OUTER_SHELL); -// upper= total_fb (&wmain[xxxplasma->nwind], t, f1+(n*df), f2, FB_FULL, OUTER_SHELL); -// test= total_fb (&wmain[xxxplasma->nwind], t, f1+((n-1)*df), f1+((n)*df), FB_FULL, OUTER_SHELL); -// printf ("FB_TEST1 break %e lower %e upper %e total %e\n",f1+(n*df),lower,upper,lower+upper); -// printf ("FB_TEST2 %e to %e = %e\n",f1+((n-1)*df), f1+((n)*df),test); -// printf ("FB_TEST3 %e %e \n",f1+((n-1)*df),fb (xxxplasma, t, f1+((n-1)*df), nions, FB_FULL)); -// } - - -// exit(0); - + + xxxplasma->cool_dr = total_fb (&wmain[xxxplasma->nwind], t, 0, VERY_BIG, FB_REDUCED, INNER_SHELL); /* 78b - nsh adding this line in next to calculate direct ionization cooling without generating photons */ @@ -174,11 +156,11 @@ xtotal_emission (one, f1, f2) //The first term here is the fb cooling due to macro ions and the second gives //the fb cooling due to simple ions. //total_fb has been modified to exclude recombinations treated using macro atoms. - cooling = xplasma->cool_rr; //Note: This the fb_matom call makes no use of f1 or f2. They are passed for //now in case they should be used in the future. But they could //also be removed. // (SS) + cooling = xplasma->cool_rr; xplasma->lum_lines = total_bb_cooling (xplasma, t_e); cooling += xplasma->lum_lines; /* total_bb_cooling gives the total cooling rate due to bb transisions whether they @@ -314,6 +296,7 @@ wind between freqencies f1 and f2 11sep nsh 70 Modifications in incorporate DR cooling (very approximate at the moment) 12sep nsh 73 Added a counter for adiabatic luminosity (!) 13jul nsh 76 Split up adiabatic luminosity into heating and cooling. + 17aug nsh xchanges to wind cooling to reflect the way we now **************************************************************/ @@ -336,7 +319,8 @@ wind_cooling (f1, f2) if (wmain[n].vol > 0.0) { nplasma = wmain[n].nplasma; - cool += x = cooling (&plasmamain[nplasma],plasmamain[nplasma].t_e); + cool += x = cooling (&plasmamain[nplasma],plasmamain[nplasma].t_e); //1708 - changed this call - now computes cooling rather than luminosity - also we popultate a local + //array called cool, rather than the xplasma array - this was overrwting the xplasma array with incorrect data. lum_lines += plasmamain[nplasma].lum_lines; cool_rr += plasmamain[nplasma].cool_rr; lum_ff += plasmamain[nplasma].lum_ff;