program explicit_model_IV
    implicit none
    
    integer :: ii, max_iter
    double precision :: q, k, beta, alpha
    double precision :: G, G_ref, CR_PTC, T, T_cell, T_ref
    double precision :: Voc_ref, Isc_ref, Jsc_ref, APV
    double precision :: Vt, n_min, n_max, coef, n
    double precision :: Voc, Isc, I0, IL, V, I
    double precision :: Vmp, Imp, Pmp, P, MPPT
    double precision :: Rpv
    double precision :: dI, f, df, tol

    double precision, dimension(501) :: Vdata, Idata

    !external :: plot_IV_ascii

    ! Parámetros
    q = 1.602e-19
    k = 1.381e-23
    beta  = -0.0023
    alpha =  0.0006

    G = 960  !990.0
    G_ref = 1353.0
    CR_PTC = 16.92 !20

    T_cell = 59.4  !54.85
    T = T_cell + 273.15
    T_ref = 25 + 273.15

    Jsc_ref = 11.8e-3 !16.8e-3
    APV = 12.0
    
    Voc_ref = 2.63
    Isc_ref = Jsc_ref * APV

    Vt = (k*T)/q

    n_max = 2.0
    n_min = 1.0
    coef = 0.0015
    n = n_max + (n_min - n_max)*exp(-coef*G)

    Voc = Voc_ref + ((n*k*T)/q)*log(G/G_ref) + beta*(T - T_ref)
    Isc = Isc_ref * (G/G_ref) * (1 + alpha*(T - T_ref)) * CR_PTC

    I0 = Isc / (exp(Voc/(n*Vt)) - 1.0)
    IL = Isc

    Rpv = 0.3 !0.086

    tol = 1e-8
    max_iter = 40
    Pmp = 0.0

    do ii=0, 500
        V = Voc * dble(ii)/500.0
        I = Isc

        do
            f = I - IL + I0*(exp((V + I*Rpv)/(n*Vt)) - 1.0)
            df = 1.0 + I0*exp((V + I*Rpv)/(n*Vt))*(Rpv/(n*Vt))
            dI = -f/df
            I = I + dI
            if (abs(dI) < tol) exit
        end do

        if (I < 0.0) I = 0.0
        P = V * I

        Vdata(ii+1) = V
        Idata(ii+1) = I

        if (P > Pmp) then
            Pmp = P
            Vmp = V
            Imp = I
        end if
    end do

    MPPT = Pmp * 40

    print *, " "
    print *, " Voc_ref = ", Voc_ref
    print *, " Isc_ref = ", Isc_ref

    print *, " "
    print *, " Voc = ", Voc
    print *, " Isc = ", Isc
    
    print *, " MPPT encontrado:"
    print *, " Vmp = ", Vmp
    print *, " Imp = ", Imp
    print *, " Pmp = ", Pmp

    print *, " MPPT = ", MPPT
    print *, " "

    !call plot_IV_ascii(Vdata, Idata, 501)

end program explicit_model_IV