1) changed default_opacity from "Kramers" to "OPAL"

2) renamed input parameter Rhot_Mdotzero_factor -> Rfront_Mdotzero_factor (is not important variable now: its recommended value 1 is default)
3) added input parameter DIM_front_Mdot_factor (default 2.3). At the border between hod and cold zones, there is a negative
accretion rate Mdot = - DIM_front_Mdot_factor * Mdot(Rin)
4) added input parameter DIM_front_approach (default "maxFvis") - the disc equation is solved to the point where Mdot=0 (default) or
where Mdot = - DIM_front_Mdot_factor * Mdot(Rin)
5) added input parameter scatter_by_corona (yes/no)
6) added DiscEqFailException: public std::exception to exceptions.hpp to report about non-converging solution
7) added household variables R_dotM0_before_shift, maxR_Qirr_no_role to freddi_state
8) added boost::optional<vecd> Shadow : ring is shadowed (1) or not (0, default)
9) added inline void set_Mdot_outer_boundary() to change current_.Mdot_outer_boundary
10) added inline void set_R_dotM0_before_shift() to change current_.R_dotM0_before_shift
11) added inline void set_maxR_Qirr_no_role() to change current_.maxR_Qirr_no_role
12) added functions: Rfront_Rhot(double r, double z_r), obtain_Mdot_outer_boundary(), Tirr_critical(double r, int ii), and verify_disc_mass(double tau)
13) added 2nd argument to  R_cooling_front(double r, double sigma_at_r) and v_cooling_front(double r, double sigma_at_r);
14) added functions: double R_vis_struct (double r, double z2r_at_r) : estimate the radius, to which viscous process can reach in time tau
15) changed name of function: ring_state_vertical() -> check_ring_is_cold()
16) in freddi_evolution.cpp: changed name of function: Mdot_out() -> Mdot_outer_boundary()
17) in freddi_evolution.cpp: moved all conditions to find boundary of the hot disc to function check_ring_is_cold(); 1 : ring is cold
18) removed Tirr_exceed_critical()
19) in freddi_state.cpp: added vecd& FreddiState::Shadow() ; double v_visc(double r, double z2r)
20) in freddi_state.cpp: in const vecd& FreddiState::Qx() : added (1.0 - Shad[i]) to Qx
21) in v_cooling_front() added check that if (sigma<0.) {sigma=0.0;} in the Ludwig's scheme
22) in double FreddiState::Teff_plus(double r) const
minimum temperature on the hot branch: changed from Lasota et al., A&A 486, 523–528 (2008) to  Tavleev+23
23) in function  R_cooling_front : changed result:
R()[last()] - v_cooling_front(r) * args().calc->tau; -> r - v_cooling_front(r, sigma_at_r) * args().calc->tau;
24) in freddi_state.cpp: added function double FreddiState::Rfront_Rhot - but it is not used now
because Rfront_Mdotzero_factor == 1.0 is recommended
25) added double FreddiState::Tirr_critical (double r, int ii) to determine critical level Tirr which keeps the ring hot.
26) in freddi_state.cpp: changed name of function ring_state_vertical() to   FreddiState::check_ring_is_cold() :  returns 1 for hot, 0 for cold
27) in nonlinear_diffusion.cpp: added include "exceptions.hpp"
28) in nonlinear_diffusion.cpp: added vatriables iter_sol and maxiter to keep track on the number of iterations :
if DIM_front_approach is "outflow" (accretion rate is negative at Rhot), disk equation sometime does not converge
29) in ns_evolution.cpp: added std output to exit from error: "R_dead is positive and less than R_cor"
30) in output.cpp: added to Derived values freddi->args().basic->rin
31) in output.cpp: added to BasicFreddiFileOutput::diskStructureDump: ### Mdot = freddi->Mdot_in()
32) input parameter DIM_front_approach is either "outflow" of "maxFvis" (default)

cpp/include/arguments.hpp

@@ -175,15 +175,18 @@ class DiskStructureArguments {
 		vecd operator()(const vecd& h) const override;
 	};
 public:
-	constexpr static const char default_opacity[] = "Kramers";
+	constexpr static const char default_opacity[] = "OPAL";
 	constexpr static const double default_Mdotout = 0.;
 	constexpr static const char default_boundcond[] = "Teff";
 	constexpr static const double default_Thot = 0.;
 	constexpr static const double default_Tirr2Tvishot = 0.;
-	constexpr static const double default_Rhot_Mdotzero_factor = 1.;
+	constexpr static const double default_Rfront_Mdotzero_factor = 1.;
+	constexpr static const double default_DIM_front_Mdot_factor = 2.3;
         constexpr static const char default_check_state_approach[] = "before2024";
         constexpr static const char default_check_Sigma_approach[] = "simple";
 	constexpr static const char default_check_Temp_approach[] = "const";
+	constexpr static const char default_DIM_front_approach[] = "maxFvis";
+	constexpr static const char default_scatter_by_corona[] = "yes";
 	constexpr static const char default_initialcond[] = "powerF";
 	constexpr static const char default_wind[] = "no";
 public:
@@ -195,10 +198,13 @@ class DiskStructureArguments {
 	std::string boundcond;
 	double Thot;
 	double Tirr2Tvishot;
-	double Rhot_Mdotzero_factor;
+	double Rfront_Mdotzero_factor;
+	double DIM_front_Mdot_factor;
 	std::string check_state_approach;
 	std::string check_Sigma_approach;
 	std::string check_Temp_approach;
+	std::string DIM_front_approach;
+	std::string scatter_by_corona;
 	double F0;
 	double Mdisk0;
 	double Mdot0;
@@ -219,38 +225,60 @@ class DiskStructureArguments {
 			const BasicDiskBinaryArguments &bdb_args,
 			const std::string& opacity,
 			double Mdotout,
-			const std::string& boundcond, double Thot, double Tirr2Tvishot,
-			double Rhot_Mdotzero_factor, const std::string& check_state_approach, const std::string& check_Sigma_approach, const std::string& check_Temp_approach,
+			const std::string& boundcond, 
+			double Thot, 
+			double Tirr2Tvishot,
+			double Rfront_Mdotzero_factor, 
+			double DIM_front_Mdot_factor, 
+			const std::string& check_state_approach, 
+			const std::string& check_Sigma_approach, 
+			const std::string& check_Temp_approach,
+			const std::string& DIM_front_approach,
+			const std::string& scatter_by_corona,
 			const std::string& initialcond,
 			std::optional<double> F0,
-			std::optional<double> Mdisk0, std::optional<double> Mdot0,
+			std::optional<double> Mdisk0, 
+			std::optional<double> Mdot0,
 			std::optional<double> powerorder,
-			std::optional<double> gaussmu, std::optional<double> gausssigma,
-			const std::string& wind, const pard& windparams);
+			std::optional<double> gaussmu, 
+			std::optional<double> gausssigma,
+			const std::string& wind, 
+			const pard& windparams);
 	DiskStructureArguments(
 			const std::string &opacity,
 			const OpacityRelated &oprel,
 			double Mdotout,
 			const std::string &boundcond,
 			double Thot,
 			double Tirr2Tvishot,
-			double Rhot_Mdotzero_factor,
+			double Rfront_Mdotzero_factor,
+			double DIM_front_Mdot_factor,
 			const std::string &check_state_approach,
 			const std::string &check_Sigma_approach,
 			const std::string &check_Temp_approach,
+			const std::string &DIM_front_approach,
+			const std::string &scatter_by_corona,
 			const std::string &initialcond,
 			const std::shared_ptr<InitialFFunction> initial_F_function,
-			const std::string &wind, const pard &windparams):
-			opacity(opacity), oprel(oprel),
+			const std::string &wind, 
+			const pard &windparams):
+			opacity(opacity), 
+			oprel(oprel),
 			Mdotout(Mdotout),
-			boundcond(boundcond), Thot(Thot), Tirr2Tvishot(Tirr2Tvishot),
-			Rhot_Mdotzero_factor(Rhot_Mdotzero_factor),
+			boundcond(boundcond), 
+			Thot(Thot), 
+			Tirr2Tvishot(Tirr2Tvishot),
+			Rfront_Mdotzero_factor(Rfront_Mdotzero_factor),
+			DIM_front_Mdot_factor(DIM_front_Mdot_factor),
 			check_state_approach(check_state_approach), 
 			check_Sigma_approach(check_Sigma_approach),
 			check_Temp_approach(check_Temp_approach),
+			DIM_front_approach(DIM_front_approach),
+			scatter_by_corona(scatter_by_corona),
 			initialcond(initialcond),
 			initial_F_function(initial_F_function),
-			wind(wind), windparams(windparams) {}
+			wind(wind), 
+			windparams(windparams) {}
 	inline vecd initial_F(const vecd& h) const { return (*initial_F_function)(h); }
 	inline size_t initial_first(const vecd& h) const { return initial_F_function->first(h); }
 };
@@ -328,25 +356,27 @@ class CalculationArguments {
 	constexpr static const unsigned int default_Nt_for_tau = 200;
 	constexpr static const char default_gridscale[] = "log";
 	constexpr static const unsigned short default_starlod = 3;
+	constexpr static const int default_verb_level = 0;
 public:
 	double init_time;
 	double time;
 	double tau;
 	unsigned int Nx;
 	std::string gridscale;
 	unsigned short starlod = 3;
+	int verb_level;
 	double eps;
 public:
 	CalculationArguments(
 			double inittime,
 			double time, std::optional<double> tau,
-			unsigned int Nx, const std::string& gridscale, const unsigned short starlod,
-			double eps=1e-6):
+			unsigned int Nx, const std::string& gridscale, const unsigned short starlod, int verb_level,
+		        double eps=1e-6):
 			init_time(inittime),
 			time(time),
 			tau(tau ? *tau : time / default_Nt_for_tau),
 			Nx(Nx), gridscale(gridscale), starlod(starlod),
-			eps(eps) {}
+			verb_level(verb_level), eps(eps) {}
 };
 
 

cpp/include/exceptions.hpp

@@ -10,4 +10,11 @@ class RadiusCollapseException: public std::exception {
 	}
 };
 
+class DiscEqFailException: public std::exception {
+public:
+	virtual const char* what() const noexcept override {
+		return "Iterations >=2000";
+	}
+};
+
 #endif //FREDDI_EXCEPTIONS_HPP

cpp/include/freddi_state.hpp

@@ -243,12 +243,16 @@ class FreddiState {
 	public:
 		double Mdot_out;
 		double Mdot_in_prev = -INFINITY;
+		double Mdot_outer_boundary;
 		double t;
 		size_t i_t;
 		size_t first;
 		size_t last;
 		vecd F;
 		double F_in;
+		double R_dotM0_before_shift = 0.0;
+		double maxR_Qirr_no_role = 0.0;
+		//double DIM_front_Mdot_factor = 0.0;
 		explicit CurrentState(const DiskStructure& str);
 		CurrentState(const CurrentState&) = default;
 		CurrentState& operator=(const CurrentState&) = default;
@@ -261,7 +265,7 @@ class FreddiState {
 		boost::optional<double> Mdisk;
 		boost::optional<double> Lx;
 		boost::optional<double> Mdot_wind;
-		boost::optional<vecd> W, Tph, Qx, Tph_vis, Tph_X, Tirr, Kirr, Sigma, Height;
+		boost::optional<vecd> W, Tph, Qx, Tph_vis, Tph_X, Tirr, Kirr, Sigma, Height, Shadow;
 	};
 
 protected:
@@ -304,22 +308,34 @@ class FreddiState {
 // current_
 public:
 	inline double Mdot_out() const { return current_.Mdot_out; }
+	inline double Mdot_outer_boundary() const { return current_.Mdot_outer_boundary; }
 	inline double F_in() const { return current_.F_in; }
 	inline const vecd& F() const { return current_.F; }
 	inline double t() const { return current_.t; }
 	inline size_t i_t() const { return current_.i_t; };
 	inline size_t first() const { return current_.first; }
 	inline size_t last() const { return current_.last; }
 	inline double Mdot_in_prev() const { return current_.Mdot_in_prev; }
+	inline double R_dotM0_before_shift()  const { return current_.R_dotM0_before_shift; }
+	inline double maxR_Qirr_no_role()  const  { return current_.maxR_Qirr_no_role; }
+	//inline double DIM_front_Mdot_factor()  const  { return current_.DIM_front_Mdot_factor; }
 protected:
 	inline void set_Mdot_in_prev(double Mdot_in) { current_.Mdot_in_prev = Mdot_in; }
 	inline void set_Mdot_in_prev() { set_Mdot_in_prev(Mdot_in()); }
+	inline void set_Mdot_outer_boundary(double Mdot_boundary) { current_.Mdot_outer_boundary = Mdot_boundary; }
+	inline void set_R_dotM0_before_shift(double r) { current_.R_dotM0_before_shift = r; }
+	inline void set_maxR_Qirr_no_role(double r) {current_.maxR_Qirr_no_role = r;}
+	//inline void set_DIM_front_Mdot_factor(double r) { current_.DIM_front_Mdot_factor = r; }
 	virtual IrradiatedStar::sources_t star_irr_sources();
 public:
 	inline double omega_R(double r) const { return std::sqrt(GM() / (r*r*r)); }
 	inline double omega_i(size_t i) const { return omega_R(R()[i]); }
 	virtual double Mdot_in() const;
 	virtual double Lbol_disk() const;
+	virtual double Rfront_Rhot(double r, double z_r) const;
+	virtual double obtain_Mdot_outer_boundary();
+	virtual double Tirr_critical(double r, int ii);
+	virtual void verify_disc_mass(double tau);
 	double phase_opt() const;
 // wind_
 public:
@@ -397,6 +413,7 @@ class FreddiState {
 	const vecd& Tirr();
 	const vecd& Kirr();
 	const vecd& Height();
+	const vecd& Shadow();
 	double Luminosity(const vecd& T, double nu1, double nu2) const;
 	inline double magnitude(const double lambda, const double F0) {
 		return -2.5 * std::log10(I_lambda<HotRegion>(lambda) * cosiOverD2() / F0);
@@ -425,9 +442,11 @@ class FreddiState {
 	double Sigma_minus(double r) const;
 	double Sigma_plus(double r) const;
 	double Teff_plus(double r) const;
-	double R_cooling_front(double r);
-	double v_cooling_front(double r);
-	int ring_state_vertical(const int i); // check ring hot (1) or cold (0)
+	double R_cooling_front(double r, double sigma_at_r);
+	double R_vis_struct (double r, double z2r_at_r);
+	double v_cooling_front(double r, double sigma_at_r);
+	double v_visc(double r, double z2r);
+	int check_ring_is_cold(const int i); // check ring hot (1) or cold (0)
 	int Tirr_exceed_critical(const int i); // check Tirr > Tcrit (1) or not (0)
 };
 

cpp/include/nonlinear_diffusion.hpp

@@ -7,7 +7,7 @@
 #include <exception>	// std::exception
 #include <functional>	// std::function
 #include <vector>
-
+#include <iostream>
 
 typedef std::vector<double> vecd;
 

cpp/include/ns/ns_arguments.hpp

@@ -114,7 +114,8 @@ class NeutronStarDiskStructureArguments: public DiskStructureArguments {
 			const std::string& opacity,
 			double Mdotout,
 			const std::string& boundcond, double Thot, double Tirr2Tvishot,
-			double Rhot_Mdotzero_factor, const std::string& check_state_approach, const std::string& check_Sigma_approach, const std::string& check_Temp_approach,
+			double Rfront_Mdotzero_factor, double DIM_front_Mdot_factor, const std::string& check_state_approach, const std::string& check_Sigma_approach, const std::string& check_Temp_approach,
+			const std::string& DIM_front_approach, const std::string& scatter_by_corona,
 			const std::string& initialcond,
 			std::optional<double> F0,
 			std::optional<double> Mdisk0, std::optional<double> Mdot0,

cpp/include/orbit.hpp

@@ -14,7 +14,7 @@ double efficiency_of_accretion(double kerr);
 
 class BlackHoleFunctions {
 public:
-	static inline double r_kerrISCO(const double Mx, const double kerr) { return rgToCm(r_kerrISCORg(kerr), Mx); }
+	static inline double r_kerrISCO(const double Mx, const double kerr) {   return rgToCm(r_kerrISCORg(kerr), Mx); }
 };
 
 

cpp/src/arguments.cpp

@@ -24,6 +24,8 @@ constexpr const double DiskStructureArguments::default_Mdotout;
 constexpr const char DiskStructureArguments::default_initialcond[];
 constexpr const double DiskStructureArguments::default_Thot;
 constexpr const char DiskStructureArguments::default_boundcond[];
+constexpr const char DiskStructureArguments::default_scatter_by_corona[];
+constexpr const char DiskStructureArguments::default_DIM_front_approach[];
 constexpr const double DiskStructureArguments::mu;
 constexpr const char DiskStructureArguments::default_wind[];
 
@@ -32,8 +34,11 @@ DiskStructureArguments::DiskStructureArguments(
 		const std::string& opacity,
 		double Mdotout,
 		const std::string& boundcond, double Thot, double Tirr2Tvishot,
-		double Rhot_Mdotzero_factor, const std::string& check_state_approach, const std::string& check_Sigma_approach,
+		double Rfront_Mdotzero_factor, double DIM_front_Mdot_factor, const std::string& check_state_approach, 
+		const std::string& check_Sigma_approach,
 		const std::string& check_Temp_approach,
+		const std::string& DIM_front_approach,
+		const std::string& scatter_by_corona,
 		const std::string& initialcond,
 		std::optional<double> F0,
 		std::optional<double> Mdisk0, std::optional<double> Mdot0,
@@ -47,8 +52,8 @@ DiskStructureArguments::DiskStructureArguments(
 		boundcond(boundcond),
 		Thot(Thot),
 		Tirr2Tvishot(Tirr2Tvishot),
-		Rhot_Mdotzero_factor(Rhot_Mdotzero_factor),
-		check_state_approach(check_state_approach),check_Sigma_approach(check_Sigma_approach), check_Temp_approach(check_Temp_approach),
+		Rfront_Mdotzero_factor(Rfront_Mdotzero_factor), DIM_front_Mdot_factor(DIM_front_Mdot_factor), check_state_approach(check_state_approach),check_Sigma_approach(check_Sigma_approach), check_Temp_approach(check_Temp_approach),
+		DIM_front_approach(DIM_front_approach),scatter_by_corona(scatter_by_corona),
 		initialcond(initialcond),
 		initial_F_function(initializeInitialFFunction(oprel,
 													  bdb_args, initialcond,
@@ -376,3 +381,4 @@ constexpr const unsigned int CalculationArguments::default_Nx;
 constexpr const unsigned int CalculationArguments::default_Nt_for_tau;
 constexpr const char CalculationArguments::default_gridscale[];
 constexpr const unsigned short CalculationArguments::default_starlod;
+constexpr const int CalculationArguments::default_verb_level;