# -*- coding: utf-8 -*-
""" f
Created on Sun May  3 10:47:47 2020

@author: Quirijn
Inclui 
1) Cálculo da ETa - método FAO, INCOMPLETO

""" 
import math
import sys

# Funcoes criadas

# Equação de Tetens
def tet(T, A, B, C):
    '''Retorna o valor da função de Tetens 
    T em °C, Resultado em kPa
    '''

    tet1 =  A * math.exp((B * T)/(C + T))
    return tet1

# Derivada da Equação de Tetens
def dtet(T, A, B, C):
    '''Retorna o valor da derivada da função de Tetens 
    T em °C,  Resultado em kPa/°C
    '''

    dtet1 = A * math.exp(B * T / (T + C)) * (B / (C + T) - B * T / (C + T)**2)        
    return dtet1
    
# equacao de evapotranspiracao
def EVp(T, Rn, u, e_a):
    '''Retorna a evapotranspiracao potencial em mm/dia
    '''
    gamma_a = 0.067290539
    R_u = 0.8 * Rn / 1000   # o /1000 é pra transf kJ em MJ
    
    dv = dtet(T, A, B, C)
    e_sat = tet(T, A, B, C)
    
    EVpot_numerador = 0.408*dv*R_u + gamma_a*u*(e_sat-e_a)*900/(T + 273)
    EVpot_denominador = dv + gamma_a*(1 + 0.34*u)
    EVpot = EVpot_numerador/EVpot_denominador
    return EVpot

def fator_p(GrupoCultura, ETc):
    #Retorna o valor do fator p (AFD/AD) 
    
    a1 = -.00723
    a2 = -.04643
    b1 = .142589
    b2 = .453571
    a = a1 * GrupoCultura + a2
    b = b1 * GrupoCultura + b2
    
    p = a * ETc + b
    #print(p)
    return p

""" ************************************
*** Constantes 
************************************ """
# Eq de Tetens
A = .61078
B = 17.27
C = 237.3

""" ************************************
*** Abrir e ler o arquivo do cenário 
************************************ """

RunFile = "C:\Python_Curso\CropSim.run"
f = open(RunFile , 'r')
file = f.read()               # importa o arquivo inteiro numa string file
list_str = file.splitlines()  # separa o string file por linhas
SoilFile = list_str[5]
CropFile = list_str[8]

NMetFiles = int(list_str[11].split()[0])
MeteoFile = []
for i in range(NMetFiles):
    MeteoFile.append(list_str[12+i])
IniDia = int(list_str[14+NMetFiles])
#print(SoilFile, MeteoFile, CropFile, IniDia)  # para testar o input

""" ********************************************
*** Abrir e ler o arquivo do solo 
******************************************** """

Zcamada, TetaINI, TetaCC, TetaPMP, ArmazCC = [0],[''],[''],[''],['']

f = open(SoilFile , 'r')
file = f.read()               # importa o arquivo inteiro numa string file
list_str = file.splitlines()  # separa o string file por linhas

NCamadas = int(list_str[5])
ZeMax = float(list_str[8])
for i in range(1,NCamadas+1):
    Zcamada.append(float(list_str[11+i].split()[0]))
    TetaINI.append(float(list_str[11+i].split()[1]))
    TetaCC.append(float(list_str[11+i].split()[2]))
    TetaPMP.append(float(list_str[11+i].split()[3]))
    ArmazCC.append(TetaCC[i] * 10 * (Zcamada[i] - Zcamada[i-1])) 
#    print(ArmazCC[i], TetaCC[i], (Zcamada[i] - Zcamada[i-1]))
#    ArmazCC[i] = TetaCC[i] * (Zcamada[i] - Zcamada[i-1])


#print(ZeMax, TetaINI, TetaCC, TetaPMP)  # para testar o input


    
""" ********************************************
*** Abrir e ler o(s) arquivo(s) meteorológico(s) 
******************************************** """

Rad, Tmin, Tmax, Vap, Vento, Chuva, Tmed = [''],[''],[''],[''],[''],[''],['']

for j in range(NMetFiles):
    f = open(MeteoFile[j] , 'r')
    file = f.read()               # importa o arquivo inteiro numa string file
    list_str = file.splitlines()  # separa o string file por linhas

    for i in range(4,369):
        a = list_str[i].split()
        Rad.append(float(a[3]))     # ---------------------------------------------------------------------------------------------------------------------------
        Tmin.append(float(a[4]))
        Tmax.append(float(a[5]))
        Vap.append(float(a[6]))
        Vento.append(float(a[7]))
        Chuva.append(float(a[8]))
    
        Tmed.append((Tmin[i-3]+Tmax[i-3])/2)

#print(Rad[1], Vento[730])   #testar input
#print(Tmed[1], Tmin[365],Tmed[365],Tmax[730])

""" ***************************************
***  Abrir e ler o arquivo de cultura   
*************************************** """

f = open(CropFile , 'r')
file = f.read()               # importa o arquivo inteiro numa string file
list_str = file.splitlines()  # separa o string file por linhas

### Parte 1 - Tabela Kc
KcDVS, Kc, Ze = [], [], []
for i in range(7,13):
    a = list_str[i].split()
    KcDVS.append(float(a[0]))
    Kc.append(float(a[1]))
    Ze.append(float(a[2]))

### Soma térmica etc.
TSum1 = float(list_str[16].split()[0])
TSum2 = float(list_str[17].split()[0])
Tb = float(list_str[18].split()[0])
MSi = float(list_str[22].split()[0])
IAFi = float(list_str[23].split()[0])
Vida = float(list_str[27].split()[0])
SLA1 = float(list_str[28].split()[0])     ################
SLA2 = float(list_str[29].split()[0])     ################
Kext = float(list_str[32].split()[0])
RUE = float(list_str[33].split()[0])

### Parte 5 - Tabela AMAX e Famax
AMAXDVS, AMAX = [],[]

for i in range(37,43):
    a = list_str[i].split()
    AMAXDVS.append(float(a[0]))
    AMAX.append(float(a[1]))

Tfamax, FAMAX = [],[]

for i in range(46,53):
    a = list_str[i].split()
    Tfamax.append(float(a[0]))
    FAMAX.append(float(a[1]))

### Parte 6 e 7
ECf = float(list_str[56].split()[0])
ECr = float(list_str[57].split()[0])
ECc = float(list_str[58].split()[0])
ECs = float(list_str[59].split()[0])
Q10 = float(list_str[63].split()[0])
fRMf = float(list_str[64].split()[0])
fRMr = float(list_str[65].split()[0])
fRMc = float(list_str[66].split()[0])
fRMs = float(list_str[67].split()[0])

### Parte 8 - Tabela Partição
FDVS, Fr, Ff, Fc, Fs = [],[],[],[],[] 

for i in range(73,80):
    a = list_str[i].split()
    FDVS.append(float(a[0]))
    Fr.append(float(a[1]))
    Ff.append(float(a[2]))
    Fc.append(float(a[3]))
    Fs.append(float(a[4]))

### checar se fatores de partição somam 1                           
for i in range(7):
    if Fr[i]+Ff[i]+Fc[i]+Fs[i] != 1:
        print('Soma da partição do DVS ', FDVS[i], 'não é igual a 1' )
        sys.exit()

### Parte 9
DVSmort = float(list_str[83].split()[0])                        
Fmort = float(list_str[84].split()[0])                           
DVSmortF = float(list_str[85].split()[0])                         
FmortF = float(list_str[86].split()[0])                            

### Parte 10                                             $$$$$$$$$$$$$$$$$$$$$$$$$
GrupoCultura = int(list_str[90].split()[0])

""" ***************************************
***                                     ***
***         CICLO PRINCIPAL             ***
***                                     ***
*** *********************************** """

DVS = 0
Dia = IniDia - 1
NDia = 0
TSum = 0
IAFmax = 0
IAF = IAFi
Massa_r = MSi * Fr[0]
Massa_f = MSi * Ff[0]
Massa_c = MSi * Fc[0]
Massa_s = MSi * Fs[0]
Teta = TetaINI                          # ****************************
DeltaIAF = ['']
FracZ, Armaz = [],[]
for n in range(NCamadas+1):
    FracZ.append(float(0.))
    Armaz.append(float(0.))
    
ChuvaTot = 0
TransTot = 0
EvapoTot = 0
ExcedTot = 0
    
    
while DVS <= 2:
    
    Dia = Dia + 1    # Conta o dia juliano
    NDia = NDia + 1  # Conta a partir de 1
    TSum = TSum + Tmed[Dia] - Tb
    # Determinar DVS
    if TSum <= TSum1:
        DVS = TSum / TSum1
    else:
        DVS = 1 + (TSum-TSum1) / TSum2 
        
    # Interpolar o SLA entre SLA1 e SLA2 
    if DVS < 1:
        SLA = SLA1 + (SLA2-SLA1)*DVS
    else:
        SLA = SLA2
    
    # Interpolar Kc e Ze para o DVS específico 
    for i in range(1,6):
        
       if DVS < KcDVS[i]:
           Zedia = (Ze[i-1]*(KcDVS[i]-DVS) + Ze[i]*(DVS-KcDVS[i-1]))/(KcDVS[i]-KcDVS[i-1])
           Kcdia = (Kc[i-1]*(KcDVS[i]-DVS) + Kc[i]*(DVS-KcDVS[i-1]))/(KcDVS[i]-KcDVS[i-1])
           if Zedia > ZeMax:
               Zedia = ZeMax
           #print(Dia, NDia, round(DVS,2), round(Zedia,3),round(Kcdia,3))
           break

    # Cálculo da ETp - ETc
    #ETp = 3  # mm/dia - valor provisório
    
    Rn = Rad[Dia]       # em kJ/m2/d             -----------------------------------------------------------------------------------------------------------
    u = Vento[Dia]            # em m/s
    e_a = Vap[Dia]          # em kPa
    T = Tmed[Dia]          # em C
    
    ETp = EVp(T, Rn, u, e_a)
    print(f'Dia: {Dia} --- EvapTransp: {round(ETp,3):0<5}')
    
    ETc = ETp * Kcdia
    Evapo = ETc * math.exp(-IAF)
    Trans = ETc - Evapo    
    
    # Fator de fracionamento da Trans entre as camadas do solo
    for n in range(1,NCamadas+1):
        if Zcamada[n] <= Zedia:
            FracZ[n] = (Zcamada[n] - Zcamada[n-1]) / Zedia*1 
        else:
            if Zcamada[n-1] > Zedia: 
                FracZ[n] = 0
            else: 
                FracZ[n] = (Zedia - Zcamada[n-1]) / Zedia 
    #print(Dia, NDia, FracZ)

    # Chama a função do cálculo do fator p                     $$$$$$$$$$$$$$$$$$$$$$$
    p = fator_p(GrupoCultura, ETc)

   
    # Calcular nova úmidade
###    INSERIR AQUI O CÁLCULO DA REDUÇÃO DE TRANSPIRAÇÃO      ***
    Excedente = 0
    for n in range(1,NCamadas+1):
        Armaz[n] = Teta[n] * 10 * (Zcamada[n] - Zcamada[n-1])
        Armaz[n] = Armaz[n] - FracZ[n]*Trans + Excedente
        if n == 1:
            Armaz[n] = Armaz[n] - Evapo + Chuva[Dia]
        if Armaz[n] > ArmazCC[n]:
            Excedente = Armaz[n] - ArmazCC[n]
            Armaz[n] = ArmazCC[n]
        else:
            Excedente = 0
            
        Teta[n] = Armaz[n] / (Zcamada[n] - Zcamada[n-1]) / 10
        
        ChuvaTot = ChuvaTot + Chuva[Dia]
        TransTot = TransTot + Trans
        EvapoTot = EvapoTot + Evapo
        ExcedTot = ExcedTot + Excedente
        
    #print(Dia, Teta, round(ChuvaTot,2), round(TransTot,1), round(EvapoTot,2), round(ExcedTot,2))
            
    # Interpolar os fatores de partição para o DVS específico 
    for i in range(1,7):
        
       if DVS < FDVS[i]:
           Frdia = (Fr[i-1]*(FDVS[i]-DVS) + Fr[i]*(DVS-FDVS[i-1]))/(FDVS[i]-FDVS[i-1])
           Ffdia = (Ff[i-1]*(FDVS[i]-DVS) + Ff[i]*(DVS-FDVS[i-1]))/(FDVS[i]-FDVS[i-1])
           Fcdia = (Fc[i-1]*(FDVS[i]-DVS) + Fc[i]*(DVS-FDVS[i-1]))/(FDVS[i]-FDVS[i-1])
           Fsdia = (Fs[i-1]*(FDVS[i]-DVS) + Fs[i]*(DVS-FDVS[i-1]))/(FDVS[i]-FDVS[i-1])
           #print(Dia, DVS, round(Frdia,3),round(Ffdia,3),round(Fcdia,3),round(Fsdia,3))
           break
           
   # Calcular fotossíntese bruta 
    Raddia = (1 - math.exp(-Kext*IAF)) * Rad[Dia] #kJ/m2
    FotoBr = Raddia * RUE * 1e-6 * 1e4 # kg/ha
    #print(Dia, Raddia, Rad[Dia], FotoBr)    
    
    # Calcular a respiração de manutenção 
    RMr = Massa_r * fRMr * (Q10 ** ((Tmed[Dia]-20)/10)) 
    RMf = Massa_f * fRMf * (Q10 ** ((Tmed[Dia]-20)/10)) 
    RMc = Massa_c * fRMc * (Q10 ** ((Tmed[Dia]-20)/10)) 
    RMs = Massa_s * fRMs * (Q10 ** ((Tmed[Dia]-20)/10)) 
    RM = RMr+RMf+RMc+RMs
    #print(Dia,FotoBr, RM, FotoBr-RM)    

   # Calcular a mortalidade de caules e raízes
    if DVS >= DVSmort:
        Massa_c = Massa_c * (1-Fmort)
        Massa_r = Massa_r * (1-Fmort)

    if DVS >= DVSmortF:
        Massa_f = Massa_f * (1-FmortF)
    
   # Calcular fotossíntese líquida e particionar
    FotoLiq = FotoBr - RM
    
    Massa_r = Massa_r + FotoLiq * Frdia * ECr
    Massa_f = Massa_f + FotoLiq * Ffdia * ECf
    Massa_c = Massa_c + FotoLiq * Fcdia * ECc
    Massa_s = Massa_s + FotoLiq * Fsdia * ECs
    
    IAF = Massa_f * SLA
    if IAF > IAFmax:
        IAFmax = IAF
    
print('IAFmax ',round(IAFmax,2), '   Massa final grãos ', round(Massa_s,1))