program PPV
  implicit none
  ! ================== Parámetros de ejemplo ==================
  double precision :: G, Tc_C, Tref_C, Itot_ref
  double precision :: Isc_ref, Voc_ref, Imp_ref, Vmp_ref
  double precision :: muIsc_ref, muVoc_ref, eg_eV
  double precision :: sercell, area_m2
  double precision :: Rs, gam, IL, Io, Voc, Isc
  double precision :: imp, vmp, pmp, ff, eta, P_ele

  ! Datos de catálogo (ejemplo)
  Isc_ref   = 6.95d0     ! [A]  corriente de cortocircuito a ref
  Voc_ref   = 3.17d0     ! [V]  voltaje de circuito abierto a ref
  Imp_ref   = 6.76d0     ! [A]  corriente en MPP a ref
  Vmp_ref   = 2.94d0     ! [V]  voltaje en MPP a ref
  Tref_C    = 25.0d0     ! [°C] temperatura de referencia (STC)
  Itot_ref  = 500000.0d0 ! [W/m2] irradiancia de referencia
  muIsc_ref =  0.0006d0  ! [A/K] coef. temp. de Isc (ejemplo)
  muVoc_ref = -0.0023d0  ! [V/K] coef. temp. de Voc (ejemplo)
  eg_eV     = 0.67d0     ! [eV]  energía de banda prohibida
  sercell   = 1.0d0      ! [-]   celdas en serie en el módulo
  area_m2   = 1.0d0      ! [m2]  área del módulo (ejemplo)

  ! Condiciones de operación
  G    = 500000.d0 ! [W/m2] irradiancia sobre el módulo
  Tc_C = 25.d0     ! [°C]   temperatura de célula (si no la sabes, usa NOCT o balance aparte)

  ! ================== Cálculo principal ==================
  call Params(G,Tc_C,Tref_C,Itot_ref, &
       Isc_ref,Voc_ref,Imp_ref,Vmp_ref,muIsc_ref,muVoc_ref,eg_eV,sercell,area_m2, &
       Rs,gam,IL,Io,Voc,Isc,imp,vmp,pmp,ff,eta)


  P_ele = 2 * pmp

  ! ================== Salida de resultados ==================
  write(*,'(a)')   '========== MPPT (modelo FV 4-parámetros) =========='
  write(*,'(a,f8.2)') 'Irradiancia [W/m2] : ', G
  write(*,'(a,f8.2)') 'Tc [°C]            : ', Tc_C
  write(*,'(a,f8.3)') 'Rs [ohm]           : ', Rs
  write(*,'(a,f8.3)') 'gamma [-]          : ', gam
  write(*,'(a,f8.3)') 'IL [A]             : ', IL
  write(*,'(a,1pe12.4)') 'Io [A]          : ', Io
  write(*,'(a,f8.3)') 'Voc [V]            : ', Voc
  write(*,'(a,f8.3)') 'Isc [A]            : ', Isc
  write(*,'(a,f8.3)') 'Vmp [V]            : ', vmp
  write(*,'(a,f8.3)') 'Imp [A]            : ', imp
  write(*,'(a,f8.3)') 'Pmp [W]            : ', pmp
  write(*,'(a,f8.3)') 'FF  [-]            : ', ff
  write(*,'(a,f8.3)') 'Eta [-]            : ', eta
  write(*,'(a,f8.3)') 'P_ele [-]          : ', P_ele

end program PPV

! ==============================================================
! Subrutina principal: calcula parámetros y MPPT
! ==============================================================
subroutine Params(G,Tc_C,Tref_C,Itot_ref, &
     Isc_ref,Voc_ref,Imp_ref,Vmp_ref,muIsc_ref,muVoc_ref,eg_eV,sercell,area_m2, &
     Rs,gam,IL,Io,Voc,Isc,imp,vmp,pmp,ff,eta)

  implicit none
  double precision, intent(in)  :: G,Tc_C,Tref_C,Itot_ref
  double precision, intent(in)  :: Isc_ref,Voc_ref,Imp_ref,Vmp_ref
  double precision, intent(in)  :: muIsc_ref,muVoc_ref,eg_eV,sercell,area_m2
  double precision, intent(out) :: Rs,gam,IL,Io,Voc,Isc,imp,vmp,pmp,ff,eta

  double precision :: Tc,Tref,q_bz,Io_ref,IL_ref,a
  integer :: ierr

  Tc    = Tc_C   + 273.15d0
  Tref  = Tref_C + 273.15d0
  q_bz  = 11604.45d0    ! Constante q/k en K/V

  ! --- Rs por bisección
  call RRS(Rs,q_bz,sercell,Tref,Imp_ref,Vmp_ref,Isc_ref,Voc_ref,muVoc_ref,muIsc_ref,eg_eV,ierr)

  ! --- Parámetros base
  gam    = q_bz*(Vmp_ref - Voc_ref + Imp_ref*Rs)/(Tref*dlog(1.d0 - Imp_ref/Isc_ref))
  a      = gam/sercell
  IL_ref = Isc_ref
  Io_ref = IL_ref/dexp(q_bz*Voc_ref/(gam*Tref))

  ! --- Escalamiento
  IL = (G/Itot_ref)*(IL_ref + muIsc_ref*(Tc - Tref))
  if (IL < 0.d0) IL = 0.d0
  Io = Io_ref*((Tc/Tref)**3)*dexp((q_bz*eg_eV/a)*((1.d0/Tref)-(1.d0/Tc)))

  if (IL > 0.d0) then
     Voc = gam*Tc/q_bz * dlog(IL/Io + 1.d0)
  else
     Voc = 0.d0
  end if
  Isc = IL

  ! --- MPPT
  call MMPPT(IL,Io,gam,Tc,q_bz,Rs,Imp_ref,muIsc_ref,Itot_ref,G,Tref,imp,vmp,pmp)

  ! --- Métricas
  if (Voc > 0.d0 .and. Isc > 0.d0) then
     ff = (vmp*imp)/(Voc*Isc)
  else
     ff = 0.d0
  end if

  if (G > 1.d-9 .and. area_m2 > 0.d0) then
     eta = pmp*12                    ! eta = pmp/(G*area_m2)
  else
     eta = 0.d0
  end if

end subroutine Params

! ==============================================================
! MPPT por Newton-Raphson
! ==============================================================
subroutine MMPPT(IL,Io,gam,T,q_bz,Rs,Imp_ref,muIsc_ref,Itot_ref,G,Tref,imp,vmp,pmp)
  implicit none
  double precision, intent(in)  :: IL,Io,gam,T,q_bz,Rs
  double precision, intent(in)  :: Imp_ref,muIsc_ref,Itot_ref,G,Tref
  double precision, intent(out) :: imp,vmp,pmp
  double precision :: imxo, imxn, f1, f1p
  integer :: it

  if (G > 0.d0) then
     imxn = (G/Itot_ref)*(Imp_ref + muIsc_ref*(T - Tref))
  else
     imp = 0.d0; vmp = 0.d0; pmp = 0.d0
     return
  end if

  do it = 1,200
     imxo = imxn
     if (imxo > IL) imxo = 0.99d0*IL
     if (Io > 0.d0) then
        f1  = imxo + (imxo - IL - Io) * ( dlog((IL - imxo + Io)/Io) - imxo*q_bz*Rs/(gam*T) ) &
              / ( 1.d0 + (IL - imxo + Io)*(q_bz*Rs/(gam*T)) )
        f1p = 2.d0 + ( dlog((IL - imxo + Io)/Io) - (q_bz*Rs*imxo/(gam*T)) ) &
              / ( (1.d0 + (IL - imxo + Io)*(q_bz*Rs/(gam*T)))**2 )
     else
        f1  = 0.d0
        f1p = 1.d0
     end if
     imxn = imxo - f1/f1p
     if (dabs(imxn - imxo) <= 1.d-3) exit
  end do

  imp = max(0.d0, min(imxn, IL))
  if (Io > 0.d0) then
     vmp = dlog(1.d0 + (IL - imp)/Io)*(T*gam/q_bz) - imp*Rs
  else
     vmp = 0.d0
  end if
  pmp = imp*vmp
  !P_ele = pmp * 1200
end subroutine MMPPT

! ==============================================================
! Cálculo de Rs por bisección (Type103)
! ==============================================================
subroutine RRS(RSeries,q_bz,sercell,Tcell_ref,Imppt_ref,Vmppt_ref,Isc_ref,Voc_ref,muVoc_ref,muIsc_ref,eg,ierr)
  implicit none
  double precision, intent(out) :: RSeries
  double precision, intent(in)  :: q_bz,sercell,Tcell_ref
  double precision, intent(in)  :: Imppt_ref,Vmppt_ref,Isc_ref,Voc_ref
  double precision, intent(in)  :: muVoc_ref,muIsc_ref,eg
  double precision :: aup, alow, anew, gamup, gamlow, gamnew
  double precision :: rsup, rslow, rsnew, ioup, iolow, ionew
  double precision :: fup, flow, fnew
  integer :: it, ierr

  ierr = 0
  rsup = ((sercell*Tcell_ref*dlog(1.d0 - Imppt_ref/Isc_ref)/q_bz)+Voc_ref - Vmppt_ref)/Imppt_ref
  aup  = 1.d0
  gamup = sercell
  ioup = Isc_ref*dexp(-q_bz*Voc_ref/(gamup*Tcell_ref))

  rslow = 0.d0
  gamlow = q_bz*(Vmppt_ref - Voc_ref)/(Tcell_ref*dlog(1.d0 - Imppt_ref/Isc_ref))
  alow   = gamlow/sercell
  iolow  = Isc_ref*dexp(-q_bz*Voc_ref/(gamlow*Tcell_ref))

  do it = 1,10000
     if (dabs(rsup - rslow) <= 5.d-4) exit
     rsnew = 0.5d0*(rsup + rslow)
     gamnew = q_bz*(Vmppt_ref - Voc_ref + Imppt_ref*rsnew)/(Tcell_ref*dlog(1.d0 - Imppt_ref/Isc_ref))
     anew = gamnew/sercell
     ionew = Isc_ref*dexp(-q_bz*Voc_ref/(gamnew*Tcell_ref))

     fup  = -muVoc_ref + (gamup/q_bz)*( dlog(1.d0 + Isc_ref/ioup) + (Tcell_ref/(Isc_ref + ioup)) &
          *(muIsc_ref - Isc_ref*((q_bz*eg/(aup*Tcell_ref*Tcell_ref)) + 3.d0/Tcell_ref)) )
     flow = -muVoc_ref + (gamlow/q_bz)*( dlog(1.d0 + Isc_ref/iolow) + (Tcell_ref/(Isc_ref + iolow)) &
          *(muIsc_ref - Isc_ref*((q_bz*eg/(alow*Tcell_ref*Tcell_ref)) + 3.d0/Tcell_ref)) )
     fnew = -muVoc_ref + (gamnew/q_bz)*( dlog(1.d0 + Isc_ref/ionew) + (Tcell_ref/(Isc_ref + ionew)) &
          *(muIsc_ref - Isc_ref*((q_bz*eg/(anew*Tcell_ref*Tcell_ref)) + 3.d0/Tcell_ref)) )

     if (flow*fnew < 0.d0) then
        rsup = rsnew
        gamup = q_bz*(Vmppt_ref - Voc_ref + Imppt_ref*rsup)/(Tcell_ref*dlog(1.d0 - Imppt_ref/Isc_ref))
        aup   = gamup/sercell
        ioup  = Isc_ref*dexp(-q_bz*Voc_ref/(gamup*Tcell_ref))
     else
        rslow = rsnew
        gamlow = q_bz*(Vmppt_ref - Voc_ref + Imppt_ref*rslow)/(Tcell_ref*dlog(1.d0 - Imppt_ref/Isc_ref))
        alow   = gamlow/sercell
        iolow  = Isc_ref*dexp(-q_bz*Voc_ref/(gamlow*Tcell_ref))
     end if
  end do

  RSeries = rsnew
  if (RSeries < 1.d-3) ierr = 1
end subroutine RRS