program solve_rs_rsh
  implicit none
  integer, parameter :: dp = kind(1.0)
  integer :: iter, maxiter
  real(dp) :: tol, alpha
  real(dp) :: Rs, Rsh    ! incógnitas
  real(dp) :: IL, I0, Im, Vm, Voc, Isc, n, kB, q, T
  real(dp) :: C, Y
  real(dp) :: f1, f2
  real(dp) :: J11, J12, J21, J22
  real(dp) :: dRs, dRsh
  real(dp) :: det

  ! --- Parámetros (ejemplo: reemplaza por los tuyos) ---
  n = 1.3
  kB = 1.38064852e-23
  q = 1.602176634e-19
  T = 32.24 + 273.15     ! K (ejemplo)
  
  Voc = 2.630            ! V (ejemplo)
  Isc = 0.168            ! A (ejemplo)
  Vm =  2.290            ! V (ejemplo)
  Im =  0.160            ! A (ejemplo)

  I0 = 3.1066e-31   !4.9e-32         ! A (ejemplo)
  IL = 3.3258       !2.879           ! A (ejemplo)

  C = n*kB*T/(q*Im)

  ! --- inicialización ---
  Rs = 0.05                  ! conjetura inicial
  Rsh = Voc / Isc               ! conjetura inicial razonable
  tol = 1.0e-8
  maxiter = 60
  alpha = 1.0

  do iter = 1, maxiter
    ! Evitar división por cero
    if (Rsh <= 0.0) then
      print *, 'Rsh no positiva, ajustando conjetura'
      Rsh = max(1.0e-6, Rsh)
    end if

    ! Evaluar Y (debe ser >0)
    Y = IL + I0 - Im - (Vm + Im*Rs)/Rsh
    if (Y <= 0.0) then
      ! aplicar damping y/o cambiar conjetura si es necesario
      alpha = alpha * 0.5
      Rs = Rs * 0.5
      Rsh = Rsh * 0.5
      if (alpha < 1.0e-6) then
        print *, 'Fallo: argumento log negativo y damping muy pequeño'
        stop
      end if
      cycle
    end if

    ! f1, f2
    f1 = Rs - C*log(Y / I0) + Vm/Im
    f2 = Rsh - Voc/Isc + Rs

    ! Jacobiano
    J11 = 1.0 + C*Im/(Rsh*Y)
    J12 = -C*(Vm + Im*Rs)/(Y*Rsh*Rsh)
    J21 = 1.0
    J22 = 1.0

    ! Resolver sistema 2x2 J * d = -F
    det = J11*J22 - J12*J21
    if (abs(det) < 1.0e-30) then
      print *, 'Determinante muy pequeño, abortando'
      stop
    end if

    dRs = (-f1*J22 - (-f2)*J12) / det
    dRsh = (J11*(-f2) - J21*(-f1)) / det

    ! Actualizar con damping
    Rs = Rs + alpha*dRs
    Rsh = Rsh + alpha*dRsh

    ! Criterio de paro
    if (max(abs(dRs), abs(dRsh)) < tol) then
      print *, 'Convergencia alcanzada en iter=', iter
      exit
    end if

    ! Si la corrección es grande o provoca Y<=0, reducir alpha
    if (Y <= 0.0 .or. abs(dRs) > 1.0e3 .or. abs(dRsh) > 1.0e6) then
      alpha = alpha * 0.5
    end if
  end do

  print '(A, F12.8)', 'Rs = ', Rs
  print '(A, F12.8)', 'Rsh = ', Rsh
end program