import java.util.*;
import java.io.*;
import static java.lang.Math.*;

/**
  * Модель кристаллизации переохлажденных капель водных растворов<br>
  * Cristallization Model (CM)<br><br>
  * Источник:<br> Чукин В.В., Платонова А.С. Кристаллизация переохлажденных капель водных растворов // Естественные и технические науки. - 2009. - №4. - С.231-236.<br><br>
  * Лицензия:<br> <a href="http://www.gnu.org/licenses/gpl.html">GNU General Public License 3</a><br><br>
  * Авторы:<br>
  * Чукин Владимир Владимирович (chukin@rshu.ru)<br>
  * Платонова Анастасия Сергеевна (n-l-s@mail.ru)<br><br>
  * Версия:<br>2009-12-10
  */


public class CM {

/** Масса молекулы воды, кг */
  public static final double mH2O = 2.99152E-26;
/** Постоянная Больцмана, Дж/(кг*К) */
  public static final double k = 1.380662E-23;
/** Постоянная Планка, Дж*с */
  public static final double h = 6.626176E-34;


/** Температура плавления льда.
  * @param aw активность воды
  * @return температура плавления льда, К
  */
  public static double meltingTemperature(double aw) {
    double T = 273.16 + 103.6*log(aw) + 15.613*pow(log(aw),2) + 54.118*pow(log(aw),3);
    if(T>0.0) return T;
    else      return 0.0;
  }


/** Равновесная активность воды.
  * @param T температура, К
  * @return равновесное значение активности воды
  */
  public static double equalWaterActivity(double T) {
    return exp((210368+131.48*T-3323730/T-41729.1*log(T))/(8.31441*T));
  }


/** Эквивалентная температура
  * @param T температура, К
  * @param aw активность воды
  * @return эквивалентная температура, К
  */
  public static double equalTemperature(double T, double aw) {
    double Aw = equalWaterActivity(T);
    double T1 = 273.16 + 103.6*log(1-aw+Aw) + 15.613*pow(log(1-aw+Aw),2) + 54.118*pow(log(1-aw+Aw),3);
    if(T1>273.15) return 273.15;
    return T1;
  }


/** Скорость гомогенного образования ледяных ядер.
  * @param T температура, К
  * @param aw активность воды
  * @return скорость гомогенного образования ледяных ядер, м<sup>-3</sup>с<sup>-1</sup>
  */
  public static double getHomogeneousIceNucleationRate(double T, double aw) {
    T = equalTemperature(T, aw);
    double rhoWater = 1000.0;
    double rhoIce = 956.58-0.144886*T;
    double sigma = 0.00025*T-0.0397875;
    double Lwi = 6.7E-21; //(-1004209.186+7766.15*T-10.5*T*T)*mH2O;
    double T0 = 273.16;
    double J0 = 2*pow(rhoWater/mH2O,2.0/3.0)*rhoWater/rhoIce*k*T/h*pow(sigma/k/T,0.5);
    double rc = 2.0*sigma*mH2O/(rhoIce*Lwi*log(T0/T));
    double w = 1.0;
    double dGmax = 4.0*PI*rc*rc*sigma*w/3.0;
    double dGact = 3.6E-20 - 7.3E-22*(T-T0);
    return J0*exp(-dGmax/k/T)*exp(-dGact/k/T);
  }


/** Температура гомогенной кристаллизации капель.
  * @param r радиус капли раствора, м
  * @param aw активность воды
  * @param dTdt скорость охлаждения, К/с
  * @return температура гомогенной кристаллизации, К
  */
  public static double homogeneousTemperature(double r, double aw, double dTdt) {
    int i;
    int N = 27316;
    double T = 273.16;
    double J_hom;
    double dT = 0.01;
    double P = 0.0;

    for(i=0;i<N;i++) {
      T = 273.16 - dT*i;
      J_hom = getHomogeneousIceNucleationRate(T, aw);
      P = P - 4.0*PI*r*r*r*J_hom*(1-P)/3.0/dTdt*dT;
      if(P>=0.5) break;
    }

    return T;
  }


/** Скорость гетерогенного образования ледяных ядер.
  * @param T температура, К
  * @param aw активность воды
  * @param alpha удельная поверхностная энергия, Дж/м<sup>2</sup>
  * @return скорость гетерогенного образования ледяных ядер, м<sup>-2</sup>с<sup>-1</sup>
  */
  public static double getHeterogeneousIceNucleationRate(double T, double aw, double alpha) {
    if(aw!=1.0) T = equalTemperature(T, aw);
    double rhoWater = 1000.0;
    double rhoIce = 956.58-0.144886*T;
    double sigma = 0.00025*T-0.0397875;
    double Lwi = 7E-23*T-1E-20; //6.7E-21; //(-1004209.186+7766.15*T-10.5*T*T)*mH2O;
    double T0 = 273.16;
    double J0 = pow(rhoWater/mH2O,2.0/3.0)*k*T/h;
//    double rc = 2.0*sigma*mH2O/(rhoIce*Lwi*log(T0/T));
    double rc = alpha*pow(mH2O/rhoWater, 2.0/3.0)/(Lwi*log(T0/T));
    double dGmax = PI*rc*alpha;
    double dGact = 3.6E-20 - 7.3E-22*(T-T0);
//    System.err.printf(Locale.US, "*** %.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\r\n", T, J0, Lwi, rc, dGmax, dGact);
    return J0*exp(-dGmax/k/T)*exp(-dGact/k/T);
  }


/** Температура гетерогенной кристаллизации капель.
  * @param r радиус капли раствора, м
  * @param aw активность воды
  * @param Na концентрация нерастворимых частиц в капле раствора, м<sup>-3</sup>
  * @param alpha удельная поверхностная энергия, Дж/м<sup>2</sup>
  * @param dTdt скорость охлаждения, К/с
  * @return температура гетерогенной кристаллизации, К
  */
  public static double heterogeneousTemperature(double r, double aw, double ra, double Na, double alpha, double dTdt) {
    int i;
    int N = 27316;
    double T = 273.16;
    double J_het;
    double dT = 0.01;
    double P = 0.0;

    for(i=0;i<N;i++) {
      T = 273.16 - dT*i;
      J_het = getHeterogeneousIceNucleationRate(T, aw, alpha);
      P = P - 16.0*PI*PI*ra*ra*Na*r*r*r*J_het*(1-P)/3.0/dTdt*dT;
//      System.out.printf(Locale.US, "%.2f\t%.2e\t%.2f\n", T, J_het, P);
      if(P>=0.5) break;
    }

    return T;
  }

}