.. _program_listing_file_src_phase.cpp:

Program Listing for File phase.cpp
==================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_phase.cpp>` (``src/phase.cpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "phase.h"
   //Conditionals for Phase Diagrams in each layer begin LINE 249
   
   PhaseDgm::PhaseDgm(string Comp_name, EOS* (*f)(double, double), int k, EOS** phase_name, double *start_pressure)
   {
     Comp_type = Comp_name;
     phase_lowP = f;
     n = k;
     if(k > 0)
       phase_list = new EOS*[k];
     if(k > 1)
       start_P = new double[k-1];  // In GPa
     
     for(int i=0; i<n; i++)
     {
       phase_list[i] = phase_name[i];
       if(i != n-1)
         start_P[i] = start_pressure[i];
     }
   }
   
   PhaseDgm::PhaseDgm(const PhaseDgm &a)
   {
     Comp_type = a.Comp_type;
     n = a.n;
   
     if(n > 0)
       phase_list = new EOS*[n];
     if(n > 1)
       start_P = new double[n-1];  // In GPa
     
     for(int i=0; i<n; i++)
     {
       phase_list[i] = a.phase_list[i];
       if(i != n-1)
         start_P[i] = a.start_P[i];
     }
   
     phase_lowP = a.phase_lowP;
   }
   
   PhaseDgm::~PhaseDgm()
   {
     if(n>0)
       delete[] phase_list;
     if(n>1)
       delete[] start_P;
   }
   
   void PhaseDgm::set_phase_highP(int k, double *start_pressure, EOS** phase_name)
   // start pressure is an array with dimension k-1.  The first phase will replace the one with the same name in the phase diagram.
   {
     if(n > 1)
       delete[] start_P;
     if(n > 0)
       delete[] phase_list;
     n = k;
     phase_list = new EOS*[k];
     if(k > 1)
       start_P = new double[k-1];
     for(int i=0; i<n; i++)
     {
       phase_list[i] = phase_name[i];
       if(i != n-1)
         start_P[i] = start_pressure[i];
     }
   }
   
   EOS* PhaseDgm::find_phase(double P, double T)
   // Pressure in microbar
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
     EOS* phase_matched;
     string first_token, s;
     if(n == 0)
       return phase_lowP(P,T);
     else if(n == 1)
     {
       phase_matched = phase_lowP(P,T);
       s = phase_matched->getEOS();
       first_token = s.substr(0,s.find(" ("));
       if(first_token == phase_list[0]->getEOS())
         return phase_list[0];
       else
         return phase_matched;
     }
     else
     {
       int i = n-2;
       while(P < start_P[i]*1E10 && i > 0)
         i--;
   
       if(P >= start_P[i]*1E10)
         return phase_list[i+1];
       else
       {
         phase_matched = phase_lowP(P,T);
         s = phase_matched->getEOS();
         first_token = s.substr(0,s.find(" ("));
         if(first_token == phase_list[0]->getEOS())
           return phase_list[0];
         else
           return phase_matched;
       }
     }
   }
   
   double checkphase(double P, void *params)
   {
     struct phase_params *p = (struct phase_params *) params;
     PhaseDgm* cmpn = p->cmpn;
     double Pu=p->Pu, Pl=p->Pl;
     
     if (P < Pl || P > Pu)
     {
       cout<<"Pressure "<<P/1E10<<" is outside the bound between "<<Pl<<" and "<<Pu<<" when solving the pressure of phase boundary for component "<<cmpn->getname()<<" from phase "<<p->Phase->getEOS()<<endl;
       return numeric_limits<double>::quiet_NaN();
     }
     
     double Tt;
     if (P == Pu)          // In case the output Tt is very close but not exact Tu or Tl.
       Tt = p->Tu;
     else if (P == Pl)
       Tt = p->Tl;
     else
       Tt = (p->Tl*(Pu-P)+p->Tu*(P-Pl))/(Pu-Pl);           // Linear interpolate temperature at P
   
     if (cmpn->find_phase(P, Tt) == p->Phase)
       return 1.;
     else
       return -1.;
   }
   
   EOS* PhaseDgm::find_phase_boundary(double Pl, double Pu, double Tl, double Tu, bool inward, double &Po, double &To, double &rhoo, double &Pn, double &Tn, double &rhon)
   // Used when integrate adiabatic profile across the phase boundary.  Given the pressure in cgs, integration direction, the lower and upper pressure and temperature limit (at previous and next integral step depends on the integration direction).  Return the pressure, temperature, density at old (o) phase boundary and new (n) phase boundary.
   {
     EOS* new_phase;
     double P1, P2;
     
     int iter = 0, max_iter = 100;
     int status;
     const gsl_root_fsolver_type *TPL = gsl_root_fsolver_bisection;
   
     gsl_root_fsolver *s = gsl_root_fsolver_alloc (TPL);
     gsl_function F;
   
     EOS* old_phase;
     if (inward)
       old_phase = find_phase (Pl, Tl);
     else
       old_phase = find_phase (Pu, Tu);
       
     struct phase_params params = {Pl, Pu, Tl, Tu, old_phase, this};
     F.function = &checkphase;
     F.params = &params;
   
     status = gsl_root_fsolver_set (s, &F, Pl, Pu);
     if (status == GSL_EINVAL)
     {
       cout<<"Error: The phase transition boundary is not between the lower pressure bound "<<Pl<<" and upper pressure bound "<<Pu<<", low temperature "<<Tl<<" and high temperature "<<Tu<<" when find the phase boundary of component "<<Comp_type<<" from phase "<<old_phase->getEOS()<<" to new phase ";
   
       if (inward)
         new_phase = find_phase (Pu, Tu);
       else
         new_phase = find_phase (Pl, Tl);
       
       gsl_root_fsolver_free (s);
       Po = numeric_limits<double>::quiet_NaN();
       To = numeric_limits<double>::quiet_NaN();
       rhoo = numeric_limits<double>::quiet_NaN();
       Pn = numeric_limits<double>::quiet_NaN();
       Tn = numeric_limits<double>::quiet_NaN();
       rhon = numeric_limits<double>::quiet_NaN();
       return NULL;
     }
     else if (status != GSL_SUCCESS)
     {
       cout<<"Error: Failed to set up the phase boundary of component "<<Comp_type<<" from phase "<<old_phase->getEOS()<<endl;
       gsl_root_fsolver_free (s);
       Po = numeric_limits<double>::quiet_NaN();
       To = numeric_limits<double>::quiet_NaN();
       rhoo = numeric_limits<double>::quiet_NaN();
       Pn = numeric_limits<double>::quiet_NaN();
       Tn = numeric_limits<double>::quiet_NaN();
       rhon = numeric_limits<double>::quiet_NaN();
       return NULL;
     }
     do
     {
       iter++;
   
       status = gsl_root_fsolver_iterate (s);
   
       P1 = gsl_root_fsolver_x_lower (s);
       P2 = gsl_root_fsolver_x_upper (s);
   
       status = gsl_root_test_interval (P1, P2, 1E-10, 1E-12);
     }
     while (status == GSL_CONTINUE && iter < max_iter);
   
     if (status == GSL_CONTINUE)
     {
       cout<<"Error: Can't find the pressure of the phase boundary from "<<find_phase (Pl, Tl) -> getEOS()<<" and "<<find_phase (Pu, Tu) -> getEOS()<<" at pressure "<<Pl/1E10<<" GPa within maximum interation "<<max_iter<<endl;
       gsl_root_fsolver_free (s);
       Po = numeric_limits<double>::quiet_NaN();
       To = numeric_limits<double>::quiet_NaN();
       rhoo = numeric_limits<double>::quiet_NaN();
       Pn = numeric_limits<double>::quiet_NaN();
       Tn = numeric_limits<double>::quiet_NaN();
       rhon = numeric_limits<double>::quiet_NaN();
       return NULL;
     }
   
     if (inward)
     {
       Pn = P2;
       Po = P1;
     }
     else
     {
       Pn = P1;
       Po = P2;
     }
   
     Tn = (Tl*(Pu-Pn)+Tu*(Pn-Pl))/(Pu-Pl);         // Linear interpolate temperature at P
     To = (Tl*(Pu-Po)+Tu*(Po-Pl))/(Pu-Pl);         // Linear interpolate temperature at P
     new_phase = find_phase (Pn, Tn);
     rhon = new_phase->density (Pn, Tn, rhoo);
     rhoo = find_phase (Po, To) -> density(Po, To, rhoo);
     gsl_root_fsolver_free (s);
     return new_phase;
   }
   
   // Phase curve function in Dunaeva et al. 2010, input P in GPa, return T
   double dunaeva_phase_curve(double P, double a, double b, double c, double d, double e)
   {
     P *= 1E4;         // convert P from GPa to bar
     if (P <= 0)
     {
       cout<<"Error: P = "<<P<<" bar. Dunaeva phase curve can't handle negative pressure."<<endl;
       return numeric_limits<double>::quiet_NaN();
     }
     return a + b*P + c*log(P) + d/P + e*sqrt(P);
   }
    
   // ========== Phase Diagrams for Core  ================
   // ---------------------------------
   // Fe Default: hcp and Liquid iron
   EOS* find_phase_Fe_default(double P, double T)
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
   
     P /= 1E10;            // convert microbar to GPa
   
     // Default Core
     if( T > 12.8*P + 2424 && T > 13.7*P + 2328)   // melting curve from Dorogokupets et al. 2017, Scientific Reports. fcc and hcp Fe melting curve.
     {  
       if(P<80 || T>10000)
         return Fe_liquid2;
       else
         return Fe_liquid;
     }
     else
       return Fe_hcp;             // use hcp Iron for all regions.
   }
   
   // ---------------------------------
   // Fe_fccbcc: Includes low pressure fcc and bcc iron
   EOS* find_phase_Fe_fccbcc(double P, double T)
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
   
     P /= 1E10;            // convert microbar to GPa
       
     if(T>575+18.7*P+0.213*pow(P,2)-0.000817*pow(P,3) && P<98.5) // Anzellini et al. 2013
     {
       if(T<1991*pow((((P-5.2)/27.39)+1),1/2.38) && P<98.5)
         if(T<-41.1*P+1120)    
           return Fe_bcc;
         else
           return Fe_fcc;
       else
         return Fe_liquid;
     }
     else
     {
       if(T<3712*pow((((P-98.5)/161.2)+1),1/1.72))
         if(T<-61.2*P+1266)     // Dorogokupets et al. 2017
           return Fe_bcc;
         else  
           return Fe_hcp;
       else
         return Fe_liquid;
     }
   }
   
   // ========== Phase Diagrams for Mantle  ================
   // ---------------------------------
   // Si Default: Upper Mantle: Fo, Wds, Rwd, and liquid ; Lower Mantle: Brg, PPv
   EOS* find_phase_Si_default(double P, double T)
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
   
     P /= 1E10;            // convert microbar to GPa
     if(P > 112.5 + 7E-3*T)      // Phase transfer curve from Ono & Oganov 2005, Earth Planet. Sci. Lett. 236, 914
       return Si_PPv_Sakai;
     else if (T > 1830*pow(1+P/4.6, 0.33)) // Melting curve from Belonoshko et al. 2005 Eq. 2
       return Si_Liquid_Wolf;
     else if (P > 24.3+(-2.12E-4*T)+(-3.49E-7*pow(T, 2))) // Dorogokupets et al. 2015
       return Si_Pv;
     else if (P > 8.69+6.05E-3*T)
       return Rwd;
     else if (P > 9.45+2.76E-3*T)
       return Wds;
     else
       return Fo;
   }
   
   // ---------------------------------
   // Si Simplified: Brg, PPv, and liquid 
   EOS* find_phase_Si_simple(double P, double T)
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
   
     P /= 1E10;            // convert microbar to GPa
     if(P > 112.5 + 7E-3*T)    // Phase transfer curve from Ono & Oganov 2005, Earth Planet. Sci. Lett. 236, 914
       return Si_PPv_Sakai;
     else if (T > 1830*pow(1+P/4.6, 0.33)) // Melting curve from Belonoshko et al. 2005 Eq. 2
       return Si_Liquid_Wolf;
     else
       return Si_Pv;
   }
   
   // ---------------------------------
   // PREM tabulated mantle
   EOS* find_phase_PREM(double P, double T)
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
   
     P /= 1E10;            // convert microbar to GPa  
     if(P > 3500)
       return Si_Seager;
     else if(P > 23.83)
       return Si_BM2fit;
     else
       return Si_PREM;
   }
   
   //-----------------------------------
   // Carbon Mantle
   EOS* find_phase_C_simple(double P, double T)
   {
     if (P <= 0 || T <= 0)
       return NULL;
     P /= 1E10;      // convert microbar to GPa
     
     if (P>=970.679+(-1.52854E-2*T)+(-5.72152E-7*pow(T,2))) // Transition from Benedict et al (2018)
       return BC8;
     else if (P>=1.949+(T+273)/400)  // Transition from Kennedy and Kennedy (1976)
       return Diam;
     else
       return Graph;
   }
   
   //-----------------------------------
   // Silicon Carbide Mantle
   EOS* find_phase_SiC(double P, double T)
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
   
     P /= 1E10;  // convert microbar to GPa
   
     double transition_pressure = 69.0 - 0.001 * (T - 300.0);
   
     if (P < transition_pressure) 
       return SiC_B3_Vinet;  // Low pressure zinc blende structure (Vinet EOS)
     else 
       return SiC_B1_Vinet;  // High pressure rock salt structure (Vinet EOS)
   }
    
   // ========== Phase Diagrams for Hydrosphere  ================
   
   // ---------------------------------
   // H2O Water/Ice/Gas boundaries 
   EOS* find_phase_water_default(double P, double T)
   // input P in cgs
   {
     P /= 1E10;            // convert microbar to GPa
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
     //Below supercritical point, Liquid/Vapor/Ice Ih
     if(P<0.022064) 
     {
       if(P<0.022064*exp((647.096/T)*(-7.85951783*(1-T/647.096)+1.84408259*pow(1-T/647.096,1.5)-11.7866497*pow(1-T/647.096,3)+22.6807411*pow(1-T/647.096,6)))||T>647.096) //Vapor Line Wagner & PruB 22
       {
         if(T<1000)  //use simplified EOS above 1000 K, IAPWS below
           return Water_Vap_IAPWS; //IAPWS-R6-95
         else
           return vdW_H2O; //Van Der Walls Gas
       }
       else if(T<-32.26706488*pow(P,3)-143.6159774*pow(P,2)-74.97759097*P+273.1683519) //Ice Ih melt line, SeaFreeze fitted polynomial
         return IceIh_SF;
       else if(T>490)
         return Water_IAPWS;
       else
         return Water_SF;
     }
     // Ice Ih, III, Water Triple Point (includes Ice II region) SeaFreeze
     else if( P < 0.207592 )       
     {
       if(T > -32.26706488*pow(P,3)-143.6159774*pow(P,2)-74.97759097*P+273.1683519) //Ice Ih melt line, SeaFreeze fitted polynomial
       {
         if(T<490)
           return Water_SF;
         else if(T>1000) 
           return vdW_H2O;
         else
           return Water_Vap_IAPWS;
       }
       else if(P < 3.986300903e-09*pow(T,3)-1.789026292e-06*pow(T,2)+0.0009677787797*T+0.02631882383) // Ih to II transition, SeaFreeze fitted polynomial
         return IceIh_SF;
       else
         return IceII_SF;
     }
     // III, V, Water Triple Point, (Includes Ice II reigon) SeaFreeze
     else if( P < 0.350109 )       
     {
       if(T > 99.49717929*pow(P,3)-158.5333434*pow(P,2)+100.0993404*P+236.281993) //Ice III Melt line, SeaFreeze fitted polynomial
       {
         if(T>1280) //Use Brown <1280 K and Mazevet >1280K, there will be a density decrease if transitioning from Brown -> Mazevet
           return Water_sc_Mazevet;
         else
           return Water_Brown;
       }
       else if(P > 3.986300903e-09*pow(T,3)-1.789026292e-06*pow(T,2)+0.0009677787797*T+0.02631882383 && T<10.6802119*pow(P,3)-60.79880497*pow(P,2)+108.5464607*P+218.0347735) //IceII region, SeaFreeze fitted polynomial
         return IceII_SF;
       else if(P < 1.153769845e-08*pow(T,3)-7.956399761e-06*pow(T,2)-0.001642755909*T+0.1140942871) //Small remaining IceIh region, SeaFreeze
         return IceIh_SF;
       else
         return IceIII_SF;
     }
     // V, VI, Water Triple Point SeaFreeze
     else if(P < 0.634399)     
     {
       if(T > 47.24417282*pow(P,3)-120.4230091*pow(P,2)+143.9050041*P+218.5209061) //Ice V Melt line, SeaFreeze fitted polynomial
       {
         if(T>1280)
           return Water_sc_Mazevet; //there will be a density decrease if transitioning from Brown -> Mazevet 
         else
           return Water_Brown;
       }
       else if(T<8.754511801*pow(P,3)-42.36539788*pow(P,2)-114.2410848*P+294.9938885 && T<10.6802119*pow(P,3)-60.79880497*pow(P,2)+108.5464607*P+218.0347735) //IceII region, SeaFreeze fitted polynomial
         return IceII_SF;
       else if(P< 2.669414392e-08*pow(T,3)-1.786325179e-05*pow(T,2)+0.003113786951*T+0.2759414249) //Small remaining IceIII, SeaFreeze 
         return IceIII_SF;
       else
         return IceV_SF;
     }
     // Water, Ice VI, Ice VII Triple Point, SeaFreeze
     else if(P < 2.216)        
     {
       if(T>5.711742978*pow(P,3)-39.67340784*pow(P,2)+125.6893825*P+208.5585936 || T>355) //Ice VI Melt line, SeaFreeze fitted polynomial
       {
         if(T>1280)
           return Water_sc_Mazevet;
         else
           return Water_Brown;
       }
       else if(P< -3.607661344e-08*pow(T,3)+1.032010899e-05*pow(T,2)-0.002592049235*T+1.070211316 && T<8.754511801*pow(P,3)-42.36539788*pow(P,2)-114.2410848*P+294.9938885) //Remaining IceII region
         return IceII_SF;
       else if(P< -1.970275298e-08*pow(T,3)-1.858690183e-06*pow(T,2)+0.003737731993*T+0.1540994852) //Remaining IceV region
         return IceV_SF;
       else if(T>-1.4699e5+6.10791e-6*P*1e9+8.1529e3*log(P*1e9)-8.8439e-1*sqrt(P*1e9)) // IceVI-VII transition AQUA
         return IceVI_SF;
       else
         return IceVII_Bezacier;
     }
     // Ice VII, X, Water Triple Point. Region of possible transitional Ice VII' reported in Grande not used in default
     else if(P < 30.9)     
     {
       if(P<2.17+1.253*(pow(T/355,3.0)-1) || T>1023)   // Ice VII Melting curve Datchi et al. 2000 Phys. Rev. B 61, 6535
       {
         if(T>1280)
           return Water_sc_Mazevet;
         else
           return Water_Brown;
       }
       else
       {
         if(T<700) //Avoid extrpolating Bezacier to high temperatures where dT/dP_S is too high
           return IceVII_Bezacier;
         else
           return IceVII_Fei;
       }
     }
     // Phase diagram becomes more uncertain over 30.9 GPa
     else if (P<700)           
     {                 // use Ice X if T<2250 and P<700 otherwise supercritical similiar to AQUA 
       if(P<pow(10.0, exp(1.7818*pow(T/1634.6, 0.2408) + 0.8310*pow(T/1634.6, -1.0) - 0.1444*pow(T/1634.6, -3.0)) - 1.0)/1e9 || T>2250)    // Ice X Melting curve AQUA
         return Water_sc_Mazevet;
       else
         return IceX;
     }
     else
       return Water_sc_Mazevet;
   }
   
   // ---------------------------------
   // H2O Water/Ice boundaries primarily from Dunaeva et al. 2010
   //Simplified phase diagram with major phases as used in Huang et al. 22
   EOS* find_phase_water_legacy(double P, double T)
   // input P in cgs
   {
     P /= 1E10;            // convert microbar to GPa
     double Tt1, Tt2;
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
     if( P < 0.208566 )        // liquid water or Ice Ih (Dunaeva et al. 2010)
     {
       Tt1 = dunaeva_phase_curve(P, 273.0159, -0.0132, -0.1577, 0, 0.1516);
       if (!gsl_finite(Tt1))
       {
         cout<<"Error: Can't find the phase of H2O at P="<<P<<" GPa and T="<<T<<" K."<<endl;
         return NULL;
       }
       if(T > Tt1)
         return Water;
       else
         return IceIh;
     }
   
     else if( P < 0.6324 )     // liquid water or Ice II, III, V
     {
       Tt1 = dunaeva_phase_curve(P, 10.277, 0.0265, 50.1624, 0.5868, -4.3288);
       Tt2 = dunaeva_phase_curve(P, 5.0321, -0.0004, 30.9482, 1.0018, 0);
       if (!gsl_finite(Tt1) || !gsl_finite(Tt2))
       {
         cout<<"Error: Can't find the phase of H2O at P="<<P<<" GPa and T="<<T<<" K."<<endl;
         return NULL;
       }
       if( (P < 0.3501 && T > Tt1) || (P >= 0.3501 && T > Tt2))
         return Water;
       else
       {
         return Ice_Dummy;
       }
     }
     
     else if(P < 2.216)        // liquid water or Ice VI, VII
     {
       Tt1 = dunaeva_phase_curve(P, 4.2804, -0.0013, 21.8756, 1.0018, 1.0785);
       Tt2 = dunaeva_phase_curve(P, -47.8507, 0, -389.006, 0.9932, 28.8539);
       if (!gsl_finite(Tt1) || !gsl_finite(Tt2))
       {
         cout<<"Error: Can't find the phase of H2O at P="<<P<<" GPa and T="<<T<<" K."<<endl;
         return NULL;
       }
       if(T > Tt1)
         return Water;
       else if(T > Tt2)
         return IceVI_Bezacier;
         
       else
         return IceVII_Bezacier;
     }
   
     else if(P < 5.10)     // liquid water or Ice VII.
     {
       if(T>1200)          // A dummy melting curve to avoid ice VII EOS extrapolated to temperature too high.
         return Water_sc_Mazevet;
       else
         return IceVII_Bezacier;
     }
     else if(P < 30.9)     // liquid water or Ice VII'.
     {
       if(T>1200)          // A dummy melting curve to avoid Ice VII EOS extrapolated to temperature too high.
         return Water_sc_Mazevet;
       else
         //return IceVIIp;         //Region of possible transitional Ice VII' reported in Grande not used in default
         return IceVII_Bezacier;
     }
     else if (P<700)               // Phase diagram becomes more uncertain over 30.9 GPa
     {                     // use Ice X if T<2250 and P<700 otherwise supercritical similiar to AQUA 
       if(T<2250)
         return IceX;
       else
         return Water_sc_Mazevet;
     }
     else
       return Water_sc_Mazevet;
   }
   
   
   // ---------------------------------
   // AQUA Haldemann et al. 2020 Tabulated Ice, Liquid, Vapor, Supercritical
   EOS* find_phase_water_tabulated(double P, double T)
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
   
     P /= 1E10;            // convert microbar to GPa
     return H2O_AQUA;
   }
   
   // ========== Phase Diagram for Atmosphere  ================
   // ---------------------------------
   // Ideal Gas: Isothermal for P<100 bar, the radiative/convective boundary (Nixon & Madhusudhan 2021). Adiabatic ideal gas for P > 100 bar
   EOS* find_phase_gas_default(double P, double T)
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
     else if (P < 1E8)
       return Gas_iso;
     else
       //return vdW_H2;
       return Gas;
   }
   
   // ---------------------------------
   // H/He Gas: Isothermal ideal for P<100 bar, P>100 bar: tabulated real gas, Chabrier & Debras 2021 Y=0.275
   EOS* find_phase_HHe_tabulated(double P, double T)
   {
     if (P <= 0 || T <= 0)
     {
       return NULL;
     }
     else if (P < 1E8)
       return Gas_iso;
     else
       return Gas_hhe;
   }
   
   
   PhaseDgm core("core", find_phase_Fe_default); //Phase Diagrams for Core
   PhaseDgm core1("core1", find_phase_Fe_fccbcc); 
   PhaseDgm mant("mantle", find_phase_Si_default); //Phase Diagrams for Mantle
   PhaseDgm mant1("mantle1", find_phase_Si_simple); 
   PhaseDgm mant2("mantle2", find_phase_PREM); 
   PhaseDgm mant3("mantle3", find_phase_C_simple); //Phase Diagram for Carbon Mantle
   PhaseDgm mant4("mantle4", find_phase_SiC); //Phase Diagram for Carbide Mantle
   PhaseDgm water("water", find_phase_water_default); //Phase Diagrams for Hydrosphere
   PhaseDgm water1("water1", find_phase_water_tabulated);
   PhaseDgm water2("water2", find_phase_water_legacy);
   PhaseDgm atm("atm", find_phase_gas_default); //Phase Diagrams for Atmosphere
   PhaseDgm atm1("atm1", find_phase_HHe_tabulated);
   
   EOS* find_phase(double m, double MC, double MM, double MW, double MG, double P, double T, bool inward)
   // given the accumulated mass (in g), P (in cgs) and T, return the corresponding phase
   {
     if (inward)
     {
       if(m < MC*ME)           // iron core
         return core.find_phase(P,T);
   
       else if(m < (MM+MC)*ME)     // silicon mantle
         return mant.find_phase(P,T);
   
       else if(m < (MM+MC+MW)*ME)      // hydrosphere
         return water.find_phase(P,T);
   
       else                // return the outermost existing layer 
       {
         if(MG > 1E-18)
           return atm.find_phase(P,T);
         else if(MW > 1E-18)
           return water.find_phase(P,T);
         else if(MM > 1E-18)
           return mant.find_phase(P,T);
         else
           return core.find_phase(P,T);
       }
     }
     else
     {
       if(m <= MC*ME)          // iron core
         return core.find_phase(P,T);
   
       else if(m <= (MM+MC)*ME)        // silicon mantle
         return mant.find_phase(P,T);
   
       else if(m <= (MM+MC+MW)*ME)     // hydrosphere
         return water.find_phase(P,T);
   
       else                // return the outermost existing layer 
       {
         if(MG > 1E-18)
           return atm.find_phase(P,T);
         else if(MW > 1E-18)
           return water.find_phase(P,T);
         else if(MM > 1E-18)
           return mant.find_phase(P,T);
         else
           return core.find_phase(P,T);
       }
     }
   }  
   
   
   EOS* find_phase(double m, vector<PhaseDgm> &Comp, vector<double> M, double P, double T, bool inward)
   {
     if (Comp.size() != M.size())
     {
       cout<<"Error. The length of input vectors doesn't match."<<endl;
       exit(1);
     }
   
     double mass_layers = 0;
     int n_comp = Comp.size();
   
     for (int i=0; i<n_comp; i++)
     {
       if (M[i]*ME<1)
         continue;
       mass_layers += M[i];
       if (inward)
       {
         if (m < mass_layers*ME)
           return Comp[i].find_phase(P, T);
       }
       else
         if (m <= mass_layers*ME)
           return Comp[i].find_phase(P, T);
     }
     // if nothing matched, return the outermost existing layer
     return Comp.back().find_phase(P, T);
   }