using System;
using System.Collections.Generic;
using System.Text;

namespace GeoFly
{
    /// <summary>
    /// using Penman-Monteith formula to calculate PET
    /// </summary>
    public class PETInPM:PET
    {
        /// <summary>
        /// 
        /// </summary>
        /// <param name="tav"></param>
        /// <param name="elev"></param>
        /// <param name="yr"></param>
        /// <param name="dn"></param>
        /// <param name="lat">默认为33</param>
        public PETInPM(double tav,double elev,int yr,int dn,double lat):base(yr,dn,lat)
        {
            dTav=tav;
	        dElev=elev;
	        dTavk = dTav + 273.15;
	        pwatvap = new WaterVapor(dTav,dElev);
        }
        
        /// <summary>
        /// 空气动力学阻抗
        /// </summary>
        /// <param name="wspeed"></param>
        /// <returns></returns>
        public double AeroDynResistance(double wspeed)
        {
            double ra=-1;
	        if(wspeed<=0) 
		        wspeed = 0.01;
	        ra = 114.0/wspeed;
	        return ra;
        }
        /// <summary>
        /// 冠层阻抗
        /// </summary>
        /// <param name="CO2">CO2 含量</param>
        /// <returns></returns>
	public double CanopyResistance(double CO2)
    {
        double result=-1;
	    result=49/(1.4-0.4*CO2/330);
	    return result;
    }
        /// <summary>
        /// 同上，只是默认CO2=330
        /// </summary>
        /// <returns></returns>
    public double CanopyResistance()
    {
        return CanopyResistance(330);
    }

	/// <summary>
	/// 净长波辐射（slrp -- 日照百分率）
	/// </summary>
	/// <param name="slrp"></param>
	/// <param name="rhmd"></param>
	/// <returns></returns>
	public override double NetLongWaveRadRHmd(double slrp,double rhmd)
    {
        double avp = 0;
	    avp = pwatvap.ActVapPressure(rhmd);
        
	    dNetLong = -(0.9*slrp+0.1)*(0.34-0.139*Math.Sqrt(avp))*Consts.Stef_Boltz*Math.Pow(dTavk,4);
	    return dNetLong;
    }
	public override double NetLongWaveRadAvp(double slrp,double avp)
    {
        if(avp<=0) avp=0.5;	//(really avp may change from 0 to 2.5 kpa)
	dNetLong =  -(0.9*slrp+0.1)*(0.34-0.139*Math.Sqrt(avp))*Consts.Stef_Boltz*Math.Pow(dTavk,4);
	return dNetLong;
    }

	/// <summary>
	/// 净长波辐射（slrg -- 日照时数）
	/// </summary>
	/// <param name="slrg"></param>
	/// <param name="rhmd"></param>
	/// <returns></returns>
	public double NetLongWaveRadiationRHmd(double slrg,double rhmd)
    {
        double avp = 0;
	    avp = pwatvap.ActVapPressure(rhmd);
	    dNetLong = -(0.9*slrg+0.1)*(0.34-0.139*Math.Sqrt(avp))*Consts.Stef_Boltz*Math.Pow(dTavk,4);
	    return dNetLong;
    }
	public double NetLongWaveRadiationAvp(double slrg,double avp)
    {
    	if(avp<=0)
            avp=0.5;	//(really avp may change from 0 to 2.5 kpa)
	    dNetLong =  -(0.9*slrg+0.1)*(0.34-0.139*Math.Sqrt(avp))*Consts.Stef_Boltz*Math.Pow(dTavk,4);
	    return dNetLong;
    }

	/// <summary>
	/// 潜在蒸散量
	/// </summary>
	/// <param name="wspeed"></param>
	/// <param name="rhmd"></param>
	/// <param name="G">默认为0</param>
	/// <returns></returns>

	public double PETInPMByRHmd(double wspeed,double rhmd,double G)
    {
        double dPet=0;
	    double Hnet = NetRadiation();
	    double dlta = pwatvap.TmpVapCurveSlp();
	    double lhv = pwatvap.LatHeatVapor();
	    double comb = CombineConst();
	    double svp = pwatvap.SatuVapPressure();
	    double avp = pwatvap.ActVapPressure(rhmd);
	    double psy = pwatvap.PsychroConst();
	    double rc = CanopyResistance();
	    double ra = AeroDynResistance(wspeed);
	    dEnergy = (dlta*(Hnet-G))/(dlta+psy*(1+rc/ra))/lhv;
	    dAero = (psy*comb*(svp-avp)/ra)/(dlta+psy*(1+rc/ra))/lhv;
	    dPet = dEnergy + dAero;
	    return dPet;
    }
        /// <summary>
        /// 
        /// </summary>
        /// <param name="dn"></param>
        /// <param name="wspeed"></param>
        /// <param name="rhmd"></param>
        /// <param name="G">默认为0</param>
        /// <returns></returns>
	public double PETInPMByRHmd(int dn,double wspeed,double rhmd,double G)
    {
        double dPet = 0;
	    double Hnet = NetRadiation();
	    double dlta = pwatvap.TmpVapCurveSlp();
	    double lhv = pwatvap.LatHeatVapor();
	    double comb = CombineConst();
	    double svp = pwatvap.SatuVapPressure();
	    double avp = pwatvap.ActVapPressure(rhmd);
	    double psy = pwatvap.PsychroConst(dn);			//改进点，增加dn来模拟日的变化情况
	    double rc = CanopyResistance();
	    double ra = AeroDynResistance(wspeed);
	    dEnergy = (dlta*(Hnet-G))/(dlta+psy*(1+rc/ra))/lhv;
	    dAero = (psy*comb*(svp-avp)/ra)/(dlta+psy*(1+rc/ra))/lhv;
	    dPet = dEnergy + dAero;
	    return dPet;
    }
    /// <summary>
    /// 
    /// </summary>
    /// <param name="wspeed"></param>
    /// <param name="avp"></param>
    /// <param name="G">默认为0</param>
    /// <returns></returns>
	public double PETByRAVP(double wspeed,double avp,double G)
    {
        double svp = pwatvap.SatuVapPressure();
	    double rhmd = avp/svp;
	    return PETInPMByRHmd(wspeed,rhmd,G);
    }

        /// <summary>
        /// 净辐射
        /// </summary>
        /// <returns></returns>
	public double NetRadiation()
    {
        return dNetShort+dNetLong;
    }


	public double CombineConst()
    {
        double cc=-1;
	    cc = 1710-6.85*dTav;
	    return cc;
    }


	public WaterVapor pwatvap;
    }}