Permalink
Browse files

Fix zero Sigma from commit dd26085

  • Loading branch information...
hombit committed Mar 28, 2017
1 parent c5b2df9 commit ca93dc025fb762092cd18c833053ff66727c3f32
Showing with 40 additions and 37 deletions.
  1. +40 −37 freddi.cpp
@@ -464,12 +464,14 @@ int main(int ac, char *av[]){
output_sum << "# --rout hadn't been specified, tidal radius " << r_out/solar_radius << " Rsun was used" << std::endl;
}

vector<double> W(Nx, 0.);
W = wunc(h, F, 1, Nx-1);

for( int i_t = 0; i_t <= Time/tau + 0.001; ++i_t ){
const double t = i_t * tau;
// cout << t/DAY << endl;

vector<double> W(Nx, 0.), Tph(Nx, 0.), Tph_vis(Nx, 0.), Tph_X(Nx, 0.), Tirr(Nx,0.), Sigma(Nx, 0.), Height(Nx, 0.);
vector<double> Tph(Nx, 0.), Tph_vis(Nx, 0.), Tph_X(Nx, 0.), Tirr(Nx,0.), Sigma(Nx, 0.), Height(Nx, 0.);

Mdot_in_prev = Mdot_in;
Mdot_in = ( F.at(1) - F.at(0) ) / ( h.at(1) - h.at(0) );
@@ -496,42 +498,6 @@ int main(int ac, char *av[]){
}

const double Lx = Luminosity( R, Tph_X, nu_min, nu_max, 100 ) / pow(fc, 4.);

int ii = Nx;
if (bound_cond_type == "MdotOut"){
Mdot_out = - kMdot_out * Mdot_in;
do{
ii--;
} while( Sigma.at(ii) < Sigma_hot_disk(R[ii]) );

} else if (bound_cond_type == "fourSigmaCrit"){
do{
ii--;
// Equation from Menou et al. 1999. Factor 4 is from their fig 8 and connected to point where Mdot = 0.
} while( Sigma.at(ii) < 4 * Sigma_hot_disk(R[ii]) );
} else if ( bound_cond_type == "Teff" ){
do{
ii--;
} while( Tph.at(ii) < T_min_hot_disk );
} else if ( bound_cond_type == "Tirr" ){
if ( Mdot_in >= Mdot_in_prev and ( initial_cond_shape == "power" or initial_cond_shape == "sinusgauss" ) ){
do{
ii--;
} while( Tph.at(ii) < T_min_hot_disk );
} else{
do{
ii--;
} while( Tirr.at(ii) < T_min_hot_disk );
}
} else{
throw po::invalid_option_value(bound_cond_type);
}

if ( ii < Nx-1 ){
Nx = ii+1;
// F.at(Nx-2) = F.at(Nx-1) - Mdot_out / (2.*M_PI) * (h.at(Nx-1) - h.at(Nx-2));
h.resize(Nx);
}

const double mU = -2.5 * log10( I_lambda(R, Tph, lambdaU) * cosiOverD2 / irr0U );
const double mB = -2.5 * log10( I_lambda(R, Tph, lambdaB) * cosiOverD2 / irr0B );
@@ -599,6 +565,43 @@ int main(int ac, char *av[]){
cout << er.what() << endl;
break;
}


int ii = Nx;
if (bound_cond_type == "MdotOut"){
Mdot_out = - kMdot_out * Mdot_in;
do{
ii--;
} while( Sigma.at(ii) < Sigma_hot_disk(R[ii]) );

} else if (bound_cond_type == "fourSigmaCrit"){
do{
ii--;
// Equation from Menou et al. 1999. Factor 4 is from their fig 8 and connected to point where Mdot = 0.
} while( Sigma.at(ii) < 4 * Sigma_hot_disk(R[ii]) );
} else if ( bound_cond_type == "Teff" ){
do{
ii--;
} while( Tph.at(ii) < T_min_hot_disk );
} else if ( bound_cond_type == "Tirr" ){
if ( Mdot_in >= Mdot_in_prev and ( initial_cond_shape == "power" or initial_cond_shape == "sinusgauss" ) ){
do{
ii--;
} while( Tph.at(ii) < T_min_hot_disk );
} else{
do{
ii--;
} while( Tirr.at(ii) < T_min_hot_disk );
}
} else{
throw po::invalid_option_value(bound_cond_type);
}

if ( ii < Nx-1 ){
Nx = ii+1;
// F.at(Nx-2) = F.at(Nx-1) - Mdot_out / (2.*M_PI) * (h.at(Nx-1) - h.at(Nx-2));
h.resize(Nx);
}
}

delete oprel;

0 comments on commit ca93dc0

Please sign in to comment.