updated freddi_evolution.cpp to include condition that the hot disk c…

…annot exist for Sigma< Sigma_crit_hot (=Sigma+) while Tirr > Tcrit.

this condition ensures fast convergence to small disc radius in the case without irradiation
in FreddiEvolution::truncateOuterRadius()

cpp/include/freddi_state.hpp

@@ -424,6 +424,7 @@ class FreddiState {
 	double Mdot_wind();
 	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);
 };

cpp/src/freddi_evolution.cpp

@@ -43,26 +43,51 @@ void FreddiEvolution::truncateOuterRadius() {
 	}
 
 	auto ii = last() + 1;
+
 	if (Tirr().at(last()) / Tph_vis().at(last()) < args().disk->Tirr2Tvishot) {
 	// when irradiation is not important
-	// hot disc extends as far as Sigma>Sigma_max_cold(alpha_cold) and not farther than R_cooling_front and Tirr <= Thot 
+	// hot disc extends as far as Sigma>Sigma_max_cold(alpha_cold) and not farther than R_cooling_front and Tirr >= Thot  and Teff_vis >=Teff_plus (condition below is the opposite)
+	// or hot disk cannot exist for Sigma< Sigma_crit_hot (=Sigma+) while Tirr > Tcrit: this condition ensures fast
+	    // convergence to small disc radius in the case without irradiation
 		do {
 			ii--;
 			if (ii <= first()) throw RadiusCollapseException();
-		  } while( ( R().at(ii) > R_cooling_front ( R().at(ii)) ) && ( Sigma().at(ii) < Sigma_minus(R().at(ii)) ) && ( Tirr().at(ii) < args().disk->Thot ) );
+		  } while( 
+		        // conditions for cold zone:
+			(
+			    ( R().at(ii) > R_cooling_front ( R().at(ii)) ) 
+			&& ( Sigma().at(ii) < Sigma_minus(R().at(ii)) ) 
+			&& ( Tirr().at(ii) < args().disk->Thot ) 
+			) 
+			|| (
+			    ( Sigma().at(ii) < Sigma_plus(R().at(ii)) ) 
+			    && ( Tirr().at(ii) < args().disk->Thot ) 
+			)
+			// next line is added according to Tavleev+23 results: it overrides the line above if Thot=10000. In fact, the line above should be deleted.
+		//	&& ( Tirr().at(ii) < 9040. - 2216. ) 
+			// next line is added to prevent cooling wave when effective viscous temperature is too high:
+			//&& ( Tph_vis().at(ii) < Teff_plus(R().at(ii)) )
+		);
 		//} while( ( R().at(ii) > R_cooling_front ( R().at(ii)) ) && ( Sigma().at(ii) < Sigma_minus(R().at(ii)) ) );
 	} else if (args().disk->boundcond == "Teff") {
 	// irradiation is important, the boundary is at fixed Teff
 		do {
 			ii--;
 			if (ii <= first()) throw RadiusCollapseException();
 		} while( Tph().at(ii) < args().disk->Thot );
+
+
 	} else if (args().disk->boundcond == "Tirr") {
 	// irradiation is important, the boundary is at fixed Tir
 		do {
 			ii--;
 			if (ii <= first()) throw RadiusCollapseException();
-		} while( Tirr().at(ii) < args().disk->Thot );
+		//} while( Tirr().at(ii) < 9040. - 2216.* Tph_vis().at(ii) / Tirr().at(ii) ); 
+		    // according to Tavleev+23 results; it is applicable in quasi-stationary disc (without cooling wave)
+// 		    std::cout << Tirr().at(ii) << "\n" << std::endl;
+// 		    std::cout << args().disk->Thot << "\n" << std::endl;
+		// old: 
+		 } while( Tirr().at(ii) < args().disk->Thot );
 	} else{
 		throw std::invalid_argument("Wrong boundcond");
 	}

cpp/src/freddi_state.cpp

@@ -222,6 +222,7 @@ const vecd& FreddiState::Tirr() {
 		}
 		opt_str_.Tirr = std::move(x);
 	}
+
 	return *opt_str_.Tirr;
 }
 
@@ -696,6 +697,14 @@ double FreddiState::Sigma_plus(double r) const {
 		* std::pow(args().basic->Mx / GSL_CONST_CGSM_SOLAR_MASS, -0.37);
 }
 
+//@@@@@ 28/09/2023
+double FreddiState::Teff_plus(double r) const {
+	// Lasota et al., A&A 486, 523–528 (2008), Eq A.1, DOI: 10.1051/0004-6361:200809658
+	return 6890 *  std::pow(r / 1e10, -0.09) 
+		* std::pow(args().basic->Mx / GSL_CONST_CGSM_SOLAR_MASS, 0.03);
+}
+
+
 double FreddiState::v_cooling_front(double r) {
         // The cooling-front velocity depends on the ratio between the current Sigma and critical Sigmas
         // Ludwig et al., A&A 290, 473-486 (1994), section 3