Commit 129d024b authored by Julien Lesgourgues's avatar Julien Lesgourgues
Browse files

bug fixes and small mistakes corrected (curved lensing, nCl, annihilation, dens-lens correlations)

* 2.4:
  lensing potential in curved space: corrected a bug in the previous bug fix, and checked that everything OK now
  corrected the reading of the flag on the spot in input module
  Previous fix for longstanding issue fixed one bug but introduced another. This is the real fix.
  corrected something in reading annihilation parameters
  added few comments for recfast parameter definitions
  added few comments and tests to the DM annihilation part
  updated version number to 2.4.3
  Corrected long-standing rare bug in ndf15
  output all dens[i]-lens[j] correlations, not just i<j
  solved a bug: last_index_back and last_index_thermo not initialised in perturb_solve()
  updated hydrogen fraction like in Galli et al 2013
  updated heating fraction using Galli et al 2013
  corrected same exponents in annihilation part of hyrec
  corrected some flag and exponent typos in annihilation part
  Two remaining lines had same bug which could lead to seg-fault on 32 bit machines.
  corrected error in evolution function s impact on lensing
  added a comment concerning delta_m that needs to be generalised to ddm
  corrected lensing potential expression in cruved space
  corrected several minor typos affecting nCl but by negligible amount
  fixed missing curvature factor relevant only for nCl
parents 20aed25e 3d87ee31
This diff is collapsed.
This diff is collapsed.
......@@ -15,7 +15,7 @@
#ifndef __COMMON__
#define __COMMON__
#define _VERSION_ "v2.4.2"
#define _VERSION_ "v2.4.3"
#define _TRUE_ 1 /**< integer associated to true statement */
#define _FALSE_ 0 /**< integer associated to false statement */
......@@ -438,25 +438,21 @@ struct precision
double recfast_z_He_1; /**< down to which redshift Helium fully ionized */
double recfast_delta_z_He_1; /**< z range over which transition is smoothed */
double recfast_z_He_2; /**< down to which redshift first Helium recombination
not complete */
double recfast_z_He_2; /**< down to which redshift first Helium recombination not complete */
double recfast_delta_z_He_2; /**< z range over which transition is smoothed */
double recfast_z_He_3; /**< down to which redshift Helium singly ionized */
double recfast_delta_z_He_3; /**< z range over which transition is smoothed */
double recfast_x_He0_trigger; /**< below which Helium ionization fraction start using
full equation for Helium */
double recfast_x_He0_trigger; /**< value below which recfast uses the full equation for Helium */
double recfast_x_He0_trigger2; /**< a second threshold used in derivative routine */
double recfast_x_He0_trigger_delta; /**< x_He range over which transition is smoothed */
double recfast_x_H0_trigger; /**< below which Helium ionization fraction start using
full equation for Helium */
double recfast_x_H0_trigger; /**< value below which recfast uses the full equation for Hydrogen */
double recfast_x_H0_trigger2; /**< a second threshold used in derivative routine */
double recfast_x_H0_trigger_delta; /**< x_H range over which transition is smoothed */
double recfast_H_frac; /**< governs time at which full equation of evolution
for Tmat is used */
double recfast_H_frac; /**< governs time at which full equation of evolution for Tmat is used */
FileName hyrec_Alpha_inf_file;
FileName hyrec_R_inf_file;
......
......@@ -1169,29 +1169,41 @@ int input_read_parameters(
/* energy injection parameters from CDM annihilation/decay */
class_read_double("annihilation",pth->annihilation);
class_read_double("decay",pth->decay);
class_read_double("annihilation_variation",pth->annihilation_variation);
class_read_double("annihilation_z",pth->annihilation_z);
class_read_double("annihilation_zmax",pth->annihilation_zmax);
class_read_double("annihilation_zmin",pth->annihilation_zmin);
class_read_double("annihilation_f_halo",pth->annihilation_f_halo);
class_read_double("annihilation_z_halo",pth->annihilation_z_halo);
class_call(parser_read_string(pfc,
"on the spot",
&(string1),
&(flag1),
errmsg),
errmsg,
errmsg);
if (pth->annihilation > 0.) {
if ((flag1 == _TRUE_) && ((strstr(string1,"y") != NULL) || (strstr(string1,"Y") != NULL))) {
pth->has_on_the_spot = _TRUE_;
}
else {
pth->has_on_the_spot = _FALSE_;
class_read_double("annihilation_variation",pth->annihilation_variation);
class_read_double("annihilation_z",pth->annihilation_z);
class_read_double("annihilation_zmax",pth->annihilation_zmax);
class_read_double("annihilation_zmin",pth->annihilation_zmin);
class_read_double("annihilation_f_halo",pth->annihilation_f_halo);
class_read_double("annihilation_z_halo",pth->annihilation_z_halo);
class_call(parser_read_string(pfc,
"on the spot",
&(string1),
&(flag1),
errmsg),
errmsg,
errmsg);
if (flag1 == _TRUE_) {
if ((strstr(string1,"y") != NULL) || (strstr(string1,"Y") != NULL)) {
pth->has_on_the_spot = _TRUE_;
}
else {
if ((strstr(string1,"n") != NULL) || (strstr(string1,"N") != NULL)) {
pth->has_on_the_spot = _FALSE_;
}
else {
class_stop(errmsg,"incomprehensible input '%s' for the field 'on the spot'",string1);
}
}
}
}
class_read_double("decay",pth->decay);
/** (c) define which perturbations and sources should be computed, and down to which scale */
ppt->has_perturbations = _FALSE_;
......@@ -3081,8 +3093,13 @@ int input_default_precision ( struct precision * ppr ) {
ppr->transfer_neglect_late_source = 400.;
ppr->l_switch_limber=10.;
// For density Cl, we recommend not to use the Limber approximation
// at all, and hence to put here a very large number (e.g. 10000); but
// if you have wide and smooth selection functions you may wish to
// use it; then 30 might be OK
ppr->l_switch_limber_for_cl_density_over_z=30.;
ppr->selection_cut_at_sigma=5.;
ppr->selection_sampling=50;
ppr->selection_sampling_bessel=20;
......
......@@ -1163,7 +1163,7 @@ int output_background(
pop->error_message);
number_of_titles = get_number_of_titles(titles);
size_data = number_of_titles*pba->bt_size;
class_alloc(data,sizeof(double*)*size_data,pop->error_message);
class_alloc(data,sizeof(double)*size_data,pop->error_message);
class_call(background_output_data(pba,
number_of_titles,
data),
......@@ -1214,7 +1214,7 @@ int output_thermodynamics(
pop->error_message);
number_of_titles = get_number_of_titles(titles);
size_data = number_of_titles*pth->tt_size;
class_alloc(data,sizeof(double*)*size_data,pop->error_message);
class_alloc(data,sizeof(double)*size_data,pop->error_message);
class_call(thermodynamics_output_data(pba,
pth,
number_of_titles,
......@@ -1499,7 +1499,7 @@ int output_open_cl_file(
}
if (psp->has_dl == _TRUE_){
for (index_d1=0; index_d1<psp->d_size; index_d1++){
for (index_d2=index_d1; index_d2<=MIN(index_d1+psp->non_diag,psp->d_size-1); index_d2++){
for (index_d2=MAX(index_d1-psp->non_diag,0); index_d2<=MIN(index_d1+psp->non_diag,psp->d_size-1); index_d2++) {
sprintf(tmp,"dens[%d]-lens[%d]",index_d1+1,index_d2+1);
class_fprintf_columntitle(*clfile,tmp,_TRUE_,colnum);
}
......
......@@ -2284,6 +2284,8 @@ int perturb_solve(
ppaw.k = k;
ppaw.ppw = ppw;
ppaw.ppw->inter_mode = pba->inter_closeby;
ppaw.ppw->last_index_back = 0;
ppaw.ppw->last_index_thermo = 0;
/** - check whether we need to print perturbations to a file for this wavenumber */
......@@ -5122,8 +5124,16 @@ int perturb_einstein(
really want gauge-dependent results) */
if (ppt->has_source_delta_m == _TRUE_) {
//ppw->delta_m += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_m/k2;
ppw->delta_m -= 2.*ppw->pvecback[pba->index_bg_H_prime]/ppw->pvecback[pba->index_bg_H] * ppw->theta_m/k2;
ppw->delta_m += 3. *ppw->pvecback[pba->index_bg_a]*ppw->pvecback[pba->index_bg_H] * ppw->theta_m/k2;
// note: until 2.4.3 there was a typo, the factor was (-2 H'/H) instead
// of (3 aH). There is the same typo in the CLASSgal paper
// 1307.1459v1,v2,v3. It came from a confusion between (1+w_total)
// and (1+w_matter)=1 [the latter is the relevant one here].
//
// note2: at this point this gauge-invariant variable is only
// valid if all matter components are pressureless and
// stable. This relation will be generalised soon to the case
// of decaying dark matter.
}
if (ppt->has_source_theta_m == _TRUE_) {
......
......@@ -1666,7 +1666,7 @@ int spectra_indices(
if ((ppt->has_cl_number_count == _TRUE_) && (ppt->has_cl_lensing_potential == _TRUE_) && (ppt->has_scalars == _TRUE_)) {
psp->has_dl = _TRUE_;
psp->index_ct_dl=index_ct;
index_ct+=(psp->d_size*(psp->d_size+1)-(psp->d_size-psp->non_diag)*(psp->d_size-1-psp->non_diag))/2;
index_ct += psp->d_size*psp->d_size - (psp->d_size-psp->non_diag)*(psp->d_size-1-psp->non_diag);
}
else {
psp->has_dl = _FALSE_;
......@@ -1729,7 +1729,7 @@ int spectra_indices(
if (psp->has_dl == _TRUE_)
for (index_ct=psp->index_ct_dl;
index_ct<psp->index_ct_dl+(psp->d_size*(psp->d_size+1)-(psp->d_size-psp->non_diag)*(psp->d_size-1-psp->non_diag))/2;
index_ct < psp->index_ct_dl+(psp->d_size*psp->d_size - (psp->d_size-psp->non_diag)*(psp->d_size-1-psp->non_diag));
index_ct++)
psp->l_max_ct[ppt->index_md_scalars][index_ct] = ppt->l_lss_max;
......@@ -2304,11 +2304,10 @@ int spectra_compute_cl(
if (_scalars_ && (psp->has_dl == _TRUE_)) {
index_ct=0;
for (index_d1=0; index_d1<psp->d_size; index_d1++) {
for (index_d2=index_d1; index_d2<=MIN(index_d1+psp->non_diag,psp->d_size-1); index_d2++) {
for (index_d2=MAX(index_d1-psp->non_diag,0); index_d2<=MIN(index_d1+psp->non_diag,psp->d_size-1); index_d2++) {
cl_integrand[index_q*cl_integrand_num_columns+1+psp->index_ct_dl+index_ct]=
primordial_pk[index_ic1_ic2]
* 0.5*(transfer_ic1_nc[index_d1] * transfer_ic2[ptr->index_tt_lensing+index_d2] +
transfer_ic1[ptr->index_tt_lensing+index_d1] * transfer_ic2_nc[index_d2])
* transfer_ic1_nc[index_d1] * transfer_ic2[ptr->index_tt_lensing+index_d2]
* factor;
index_ct++;
}
......
......@@ -308,6 +308,11 @@ int thermodynamics_init(
pth->error_message,
"annihilation parameter cannot be negative");
class_test((pth->annihilation>1.e-4),
pth->error_message,
"annihilation parameter suspiciously large (%e, while typical bounds are in the range of 1e-7 to 1e-6)",
pth->annihilation);
class_test((pth->annihilation_variation>0),
pth->error_message,
"annihilation variation parameter must be negative (decreasing annihilation rate)");
......@@ -328,6 +333,10 @@ int thermodynamics_init(
pth->error_message,
"CDM annihilation effects require the presence of CDM!");
class_test((pth->annihilation_f_halo>0) && (pth->recombination==recfast),
pth->error_message,
"Switching on DM annihilation in halos requires using HyRec instead of RECFAST. Otherwise some values go beyond their range of validity in the RECFAST fits, and the thermodynamics module fails. Two solutions: add 'recombination = HyRec' to your input, or set 'annihilation_f_halo = 0.' (default).");
class_test((pth->annihilation_f_halo<0),
pth->error_message,
"Parameter for DM annihilation in halos cannot be negative");
......@@ -1243,7 +1252,7 @@ int thermodynamics_energy_injection(
class_call(thermodynamics_onthespot_energy_injection(ppr,pba,preco,zp,&onthespot,error_message),
error_message,
error_message);
first_integrand = factor*pow(1+z,6)/pow(1+zp,5.5)*exp(2./3.*factor*(pow(1+z,1.5)-pow(1+zp,1.5)))*onthespot;
first_integrand = factor*pow(1+z,8)/pow(1+zp,7.5)*exp(2./3.*factor*(pow(1+z,1.5)-pow(1+zp,1.5)))*onthespot; // beware: versions before 2.4.3, there were rwrong exponents: 6 and 5.5 instead of 8 and 7.5
result = 0.5*dz*first_integrand;
/* other points in trapezoidal integral */
......@@ -1253,7 +1262,7 @@ int thermodynamics_energy_injection(
class_call(thermodynamics_onthespot_energy_injection(ppr,pba,preco,zp,&onthespot,error_message),
error_message,
error_message);
integrand = factor*pow(1+z,6)/pow(1+zp,5.5)*exp(2./3.*factor*(pow(1+z,1.5)-pow(1+zp,1.5)))*onthespot;
integrand = factor*pow(1+z,8)/pow(1+zp,7.5)*exp(2./3.*factor*(pow(1+z,1.5)-pow(1+zp,1.5)))*onthespot; // beware: versions before 2.4.3, there were rwrong exponents: 6 and 5.5 instead of 8 and 7.5
result += dz*integrand;
} while (integrand/first_integrand > 0.02);
......@@ -1803,6 +1812,7 @@ int thermodynamics_reionization_sample(
double dkappadtau,dkappadtau_next;
double energy_rate;
double tau;
double chi_heat;
int last_index_back;
double relative_variation;
......@@ -2028,11 +2038,14 @@ int thermodynamics_reionization_sample(
pth->error_message,
pth->error_message);
//chi_heat = (1.+2.*preio->reionization_table[i*preio->re_size+preio->index_re_xe])/3.; // old approximation from Chen and Kamionkowski
chi_heat = MIN(0.996857*(1.-pow(1.-pow(preio->reionization_table[i*preio->re_size+preio->index_re_xe],0.300134),1.51035)),1); // coefficient as revised by Slatyer et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Slatyer et al. 2013)
dTdz=2./(1+z)*preio->reionization_table[i*preio->re_size+preio->index_re_Tb]
-2.*mu/_m_e_*4.*pvecback[pba->index_bg_rho_g]/3./pvecback[pba->index_bg_rho_b]*opacity*
(pba->T_cmb * (1.+z)-preio->reionization_table[i*preio->re_size+preio->index_re_Tb])/pvecback[pba->index_bg_H]
-2./(3.*_k_B_)*energy_rate*(1.+2.*preio->reionization_table[i*preio->re_size+preio->index_re_xe])
/(3*preco->Nnow*pow(1.+z,3))/(1.+preco->fHe+preio->reionization_table[i*preio->re_size+preio->index_re_xe])
-2./(3.*_k_B_)*energy_rate*chi_heat
/(preco->Nnow*pow(1.+z,3))/(1.+preco->fHe+preio->reionization_table[i*preio->re_size+preio->index_re_xe])
/(pvecback[pba->index_bg_H]*_c_/_Mpc_over_m_*(1.+z)); /* energy injection */
/** - increment baryon temperature */
......@@ -2171,6 +2184,7 @@ int thermodynamics_recombination_with_hyrec(
param.dlna = 8.49e-5;
param.nz = (long) floor(2+log((1.+param.zstart)/(1.+param.zend))/param.dlna);
param.annihilation = pth->annihilation;
param.has_on_the_spot = pth->has_on_the_spot;
param.decay = pth->decay;
param.annihilation_variation = pth->annihilation_variation;
param.annihilation_z = pth->annihilation_z;
......@@ -2864,6 +2878,8 @@ int thermodynamics_derivs_with_recfast(
double energy_rate;
double tau;
double chi_heat;
double chi_ion_H;
int last_index_back;
ptpaw = parameters_and_workspace;
......@@ -2980,23 +2996,48 @@ int thermodynamics_derivs_with_recfast(
timeTh=(1./(preco->CT*pow(Trad,4)))*(1.+x+preco->fHe)/x;
timeH=2./(3.*preco->H0*pow(1.+z,1.5));
/************/
/* hydrogen */
/************/
if (x_H > ppr->recfast_x_H0_trigger)
dy[0] = 0.;
else {
/* equations modified to take into account energy injection from dark matter */
if (x_H > ppr->recfast_x_H0_trigger2) {
dy[0] = (x*x_H*n*Rdown - Rup*(1.-x_H)*exp(-preco->CL/Tmat))/ (Hz*(1.+z))
-energy_rate*(1.-x)/(3*n)/(_L_H_ion_*_h_P_*_c_*Hz*(1.+z)); /* energy injection (neglect fraction going to helium) */
/* Peebles' coefficient (approximated as one when the Hydrogen
ionisation fraction is very close to one) */
if (x_H < ppr->recfast_x_H0_trigger2) {
C = (1. + K*_Lambda_*n*(1.-x_H))/(1./preco->fu+K*_Lambda_*n*(1.-x)/preco->fu +K*Rup*n*(1.-x));
}
else {
C=(1. + K*_Lambda_*n*(1.-x_H))/(1./preco->fu+K*_Lambda_*n*(1.-x)/preco->fu +K*Rup*n*(1.-x));
dy[0] = ((x*x_H*n*Rdown - Rup*(1.-x_H)*exp(-preco->CL/Tmat)) *(1. + K*_Lambda_*n*(1.-x_H))) /(Hz*(1.+z)*(1./preco->fu+K*_Lambda_*n*(1.-x)/preco->fu +K*Rup*n*(1.-x)))
-energy_rate*(1.-x)/(3*n)*(1./_L_H_ion_+(1.-C)/_L_H_alpha_)/(_h_P_*_c_*Hz*(1.+z)); /* energy injection (neglect fraction going to helium) */
C = 1.;
}
/* For DM annihilation: fraction of injected energy going into
ionisation and Lya excitation */
/* - old approximation from Chen and Kamionkowski: */
//chi_ion_H = (1.-x)/3.;
/* coefficient as revised by Slatyer et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Slatyer et al. 2013): */
if (x < 1.)
chi_ion_H = 0.369202*pow(1.-pow(x,0.463929),1.70237);
else
chi_ion_H = 0.;
/* evolution of hydrogen ionisation fraction: */
dy[0] = (x*x_H*n*Rdown - Rup*(1.-x_H)*exp(-preco->CL/Tmat)) * C / (Hz*(1.+z)) /* Peeble's equation with fudged factors */
-energy_rate*chi_ion_H/n*(1./_L_H_ion_+(1.-C)/_L_H_alpha_)/(_h_P_*_c_*Hz*(1.+z)); /* energy injection (neglect fraction going to helium) */
}
/************/
/* helium */
/************/
if (x_He < 1.e-15)
dy[1]=0.;
else {
......@@ -3036,8 +3077,12 @@ int thermodynamics_derivs_with_recfast(
}
else {
/* equations modified to take into account energy injection from dark matter */
//chi_heat = (1.+2.*preio->reionization_table[i*preio->re_size+preio->index_re_xe])/3.; // old approximation from Chen and Kamionkowski
chi_heat = MIN(0.996857*(1.-pow(1.-pow(x,0.300134),1.51035)),1.); // coefficient as revised by Galli et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Galli et al. 2013)
dy[2]= preco->CT * pow(Trad,4) * x / (1.+x+preco->fHe) * (Tmat-Trad) / (Hz*(1.+z)) + 2.*Tmat/(1.+z)
-2./(3.*_k_B_)*energy_rate*(1.+2.*x)/(3*n)/(1.+preco->fHe+x)/(Hz*(1.+z)); /* energy injection */
-2./(3.*_k_B_)*energy_rate*chi_heat/n/(1.+preco->fHe+x)/(Hz*(1.+z)); /* energy injection */
}
return _SUCCESS_;
......
......@@ -2190,7 +2190,24 @@ int transfer_sources(
rescaling=0.;
}
else {
rescaling = (tau_rec-tau)/(tau0-tau)/(tau0-tau_rec);
switch (pba->sgnK){
case 1:
rescaling = sqrt(pba->K)
*sin((tau_rec-tau)*sqrt(pba->K))
/sin((tau0-tau)*sqrt(pba->K))
/sin((tau0-tau_rec)*sqrt(pba->K));
break;
case 0:
rescaling = (tau_rec-tau)/(tau0-tau)/(tau0-tau_rec);
break;
case -1:
rescaling = sqrt(-pba->K)
*sinh((tau_rec-tau)*sqrt(-pba->K))
/sinh((tau0-tau)*sqrt(-pba->K))
/sinh((tau0-tau_rec)*sqrt(-pba->K));
break;
}
// Note: until 2.4.3 there was a bug here: the curvature effects had been ommitted.
}
/* copy from input array to output array */
......@@ -2264,6 +2281,10 @@ int transfer_sources(
ptr->error_message,
ptr->error_message);
class_test(tau0 - tau0_minus_tau[0] > ppt->tau_sampling[ppt->tau_size-1],
ptr->error_message,
"this should not happen, there was probably a rounding error, if this error occured, then this must be coded more carefully");
/* resample the source at those times */
class_call(transfer_source_resample(ppr,
pba,
......@@ -2595,8 +2616,8 @@ int transfer_sources(
if (_index_tt_in_range_(ptr->index_tt_nc_lens, ppt->selection_num, ppt->has_nc_lens)) {
rescaling +=
-(2.-5.*ptr->s_bias)/2.
rescaling -=
(2.-5.*ptr->s_bias)/2.
*(tau0_minus_tau[index_tau]-tau0_minus_tau_lensing_sources[index_tau_sources])
/tau0_minus_tau[index_tau]
/tau0_minus_tau_lensing_sources[index_tau_sources]
......@@ -2630,7 +2651,7 @@ int transfer_sources(
if ((ptr->has_nz_evo_file == _TRUE_) || (ptr->has_nz_evo_analytic == _TRUE_)) {
f_evo = 2./pvecback[pba->index_bg_H]/pvecback[pba->index_bg_a]/tau0_minus_tau[index_tau]
f_evo = 2./pvecback[pba->index_bg_H]/pvecback[pba->index_bg_a]/tau0_minus_tau_lensing_sources[index_tau_sources]
+ pvecback[pba->index_bg_H_prime]/pvecback[pba->index_bg_H]/pvecback[pba->index_bg_H]/pvecback[pba->index_bg_a];
z = pba->a_today/pvecback[pba->index_bg_a]-1.;
......@@ -2976,6 +2997,10 @@ int transfer_selection_sampling(
ptr->error_message,
ptr->error_message);
class_test(tau_size <= 0,
ptr->error_message,
"should be at least one");
/* case selection == dirac */
if (tau_min == tau_max) {
class_test(tau_size !=1,
......@@ -2986,9 +3011,10 @@ int transfer_selection_sampling(
/* for other cases (gaussian, tophat...) define new sampled values
of (tau0-tau) with even spacing */
else {
for (index_tau=0; index_tau<tau_size; index_tau++) {
for (index_tau=0; index_tau<tau_size-1; index_tau++) {
tau0_minus_tau[index_tau]=pba->conformal_age-tau_min-((double)index_tau)/((double)tau_size-1.)*(tau_max-tau_min);
}
tau0_minus_tau[tau_size-1]=pba->conformal_age-tau_max;
}
return _SUCCESS_;
......@@ -4376,7 +4402,8 @@ int transfer_radial_function(
//s2 = sqrt(1.0-3.0*K/k2);
factor = 1.0;
for (j=0; j<x_size; j++)
radial_function[x_size-1-j] = factor*d2Phi[j]*rescale_argument*rescale_argument*rescale_function[j];
radial_function[x_size-1-j] = factor*absK_over_k2*d2Phi[j]*rescale_argument*rescale_argument*rescale_function[j];
// Note: in previous line there was a missing factor absK_over_k2 until version 2.4.3. Credits Francesco Montanari.
break;
}
......@@ -4566,15 +4593,15 @@ int transfer_global_selection_read(
/* infer dlog(dN/dz)/dz from dN/dz */
ptr->nz_evo_dlog_nz[0] =
(ptr->nz_evo_nz[1]-ptr->nz_evo_nz[0])
(log(ptr->nz_evo_nz[1])-log(ptr->nz_evo_nz[0]))
/(ptr->nz_evo_z[1]-ptr->nz_evo_z[0]);
for (row=1; row<ptr->nz_evo_size-1; row++){
ptr->nz_evo_dlog_nz[row] =
(ptr->nz_evo_nz[row+1]-ptr->nz_evo_nz[row-1])
(log(ptr->nz_evo_nz[row+1])-log(ptr->nz_evo_nz[row-1]))
/(ptr->nz_evo_z[row+1]-ptr->nz_evo_z[row-1]);
}
ptr->nz_evo_dlog_nz[ptr->nz_evo_size-1] =
(ptr->nz_evo_nz[ptr->nz_evo_size-1]-ptr->nz_evo_nz[ptr->nz_evo_size-2])
(log(ptr->nz_evo_nz[ptr->nz_evo_size-1])-log(ptr->nz_evo_nz[ptr->nz_evo_size-2]))
/(ptr->nz_evo_z[ptr->nz_evo_size-1]-ptr->nz_evo_z[ptr->nz_evo_size-2]);
/* to test that the file is read:
......
......@@ -2279,6 +2279,17 @@ int array_interpolate_spline_growing_closeby(
int inf,sup,i;
double h,a,b;
/*
if (*last_index < 0) {
sprintf(errmsg,"%s(L:%d) problem with last_index =%d < 0",__func__,__LINE__,*last_index);
return _FAILURE_;
}
if (*last_index > (n_lines-1)) {
sprintf(errmsg,"%s(L:%d) problem with last_index =%d > %d",__func__,__LINE__,*last_index,n_lines-1);
return _FAILURE_;
}
*/
inf = *last_index;
class_test(inf<0 || inf>(n_lines-1),
errmsg,
......
......@@ -1401,7 +1401,7 @@ int numjac(
use the column computed with this increment.
This depends on wether we are in sparse mode or not: */
if ((jac->use_sparse)&&(jac->repeated_pattern >= jac->trust_sparse)){
for(i=Ap[j-1];i<Ap[j];i++) jac->xjac[i]=nj_ws->tmp[Ai[i]];
for(i=Ap[j-1];i<Ap[j];i++) jac->xjac[i]=nj_ws->tmp[Ai[i]+1];
}
else{
for(i=1;i<=neq;i++) dFdy[i][j]=nj_ws->tmp[i];
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment