program calc_Rs
  implicit none
  double precision :: Rs
  double precision :: q_bz, sercell, Tref, Imp_ref, Vmp_ref, Isc_ref, Voc_ref
  double precision :: muVoc_ref, muIsc_ref, eg

  ! ====== Datos de ejemplo ======
  Voc_ref   = 2.630d0     ! [V]
  Isc_ref   = 0.168d0     ! [A]
  Vmp_ref   = 2.290d0     ! [V]
  Imp_ref   = 0.160d0     ! [A]
  Tref      = 25.0d0 + 273.15d0  ! [K]
  muVoc_ref = -0.0023d0   ! [V/K]
  muIsc_ref = 0.0006d0    ! [A/K]
  eg        = 1.43d0      ! [eV]
  sercell   = 1 !60.d0       ! número de celdas en serie
  q_bz      = 11604.45d0  ! [K/V]

  ! ====== Cálculo de Rs ======
  call series4pam_103(Rs, q_bz, sercell, Tref, Imp_ref, Vmp_ref, Isc_ref, Voc_ref, muVoc_ref, muIsc_ref, eg)

  ! ====== Salida ======
  write(*,'(a,f10.4)') 'Resistencia serie Rs [ohm] = ', Rs

contains

  subroutine series4pam_103(RSeries,q_bz,sercell,Tcell_ref,Imppt_ref,Vmppt_ref,Isc_ref,Voc_ref,muVoc_ref,muIsc_ref,eg)
    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 :: rsup, rslow, rsnew
    double precision :: gamup, gamlow, gamnew
    double precision :: ioup, iolow, ionew
    double precision :: fup, flow, fnew
    integer :: it

    ! Límite superior e inferior
    rsup = ((sercell*Tcell_ref*dlog(1.d0 - Imppt_ref/Isc_ref)/q_bz)+Voc_ref - Vmppt_ref)/Imppt_ref
    rslow = 0.d0

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

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

    ! Funciones en los extremos
    fup  = -muVoc_ref + (gamup/q_bz)*( dlog(1.d0 + Isc_ref/ioup) + (Tcell_ref/(Isc_ref + ioup)) &
    *(muIsc_ref - Isc_ref*((q_bz*eg/(sercell*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/(sercell*Tcell_ref*Tcell_ref)) + 3.d0/Tcell_ref)) )

    ! Bisección
    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))
       ionew  = Isc_ref*dexp(-q_bz*Voc_ref/(gamnew*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/(sercell*Tcell_ref*Tcell_ref)) + 3.d0/Tcell_ref)) )

       if (flow*fnew < 0.d0) then
          rsup = rsnew
          gamup = gamnew
          ioup  = ionew
          fup   = fnew
       else
          rslow = rsnew
          gamlow = gamnew
          iolow  = ionew
          flow   = fnew
       end if
    end do

    RSeries = rsnew
  end subroutine series4pam_103

end program calc_Rs