program main_NR_3x3
  implicit none
  Double Precision :: Tabs, Tsub, Tout
  Double Precision :: G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC
  Double Precision :: W_PV, W_abs, L, A_PV, A_abs, A_ap, A_PTC
  Double Precision :: CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs
  Double Precision :: IAM_elec, IAM_th, eta_opt, sigma
  Double Precision :: mdot, Cp, U, A_hx, P_PV
  Double Precision :: Rcond_abs, Rcond_PV, Rcond_sub
  Double Precision :: Rconv_abs, Rconv_PV, Rconv_PTC

  ! ==== Valores fijos ====
  G_GHI=1000.0d0
  G_DNI=800.0d0
  Ta=25.0d0+273.15d0
  Tsky=25.0d0+273.15d0
  Tin=70.0d0+273.15d0
  TPV=355.45d0
  TPTC=310.44d0
  W_PV=0.12d0
  W_abs=0.06d0
  L=10.0d0
  A_PV=W_PV*L
  A_abs=W_abs*L
  A_ap=1.2d0*L
  A_PTC=3.0d0*L
  CR_PTC=A_ap/A_PV
  alpha_PV=0.97d0
  alpha_abs=0.90d0
  eps_PV=0.2d0
  eps_abs=0.2d0
  IAM_elec=0.72d0
  IAM_th=0.86d0
  eta_opt=0.83d0
  sigma=5.67d-8
  mdot=0.15d0
  Cp=4187.0d0
  U=780.1d0
  A_hx=0.9425d0
  P_PV=1723.0d0
  Rcond_abs=0.00002439d0
  Rcond_PV=0.00005d0
  Rcond_sub=0.00008472d0
  Rconv_abs=0.4848d0
  Rconv_PV=0.3428d0
  Rconv_PTC=0.04336d0

  ! ==== Llamada al método Newton–Raphson ====
  call solve_NR_3x3(Tabs, Tsub, Tout, &
       G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC, &
       A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs, &
       sigma, mdot, Cp, U, A_hx, P_PV, &
       Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV)

  ! ==== Resultados ====
  print *, "--------------------------------------"
  print *, "Resultados Newton–Raphson:"
  print *, "T_abs   [K]: ", Tabs
  print *, "T_sub   [K]: ", Tsub
  print *, "T_outHTF[K]: ", Tout
  print *, "--------------------------------------"

end program main_NR_3x3


!===============================================================
subroutine solve_NR_3x3(Tabs, Tsub, Tout, &
       G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC, &
       A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs, &
       sigma, mdot, Cp, U, A_hx, P_PV, &
       Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV)

  implicit none
  Double Precision, intent(out) :: Tabs, Tsub, Tout
  Double Precision, intent(in)  :: G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC
  Double Precision, intent(in)  :: A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs
  Double Precision, intent(in)  :: eps_PV, eps_abs, sigma
  Double Precision, intent(in)  :: mdot, Cp, U, A_hx, P_PV
  Double Precision, intent(in)  :: Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV

  Integer :: it, maxit
  Double Precision :: tol, normF
  Double Precision :: F(3), dx(3), J(3,3)

  tol   = 1.0d-8
  maxit = 100

  Tabs = 355.0d0
  Tsub = 355.0d0
  Tout = 355.0d0

  do it = 1, maxit
     call residual_and_jacobian(Tabs, Tsub, Tout, F, J, &
          G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC, &
          A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs, &
          sigma, mdot, Cp, U, A_hx, P_PV, &
          Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV)

     normF = max( max(abs(F(1)),abs(F(2))), abs(F(3)) )
     if (normF < tol) exit

     call solve_3x3(J, -F, dx)

     Tabs = Tabs + dx(1)
     Tsub = Tsub + dx(2)
     Tout = Tout + dx(3)
  end do

  print *, "Iteraciones: ", it-1
  print *, "||F||_inf  : ", normF

end subroutine solve_NR_3x3


!===============================================================
subroutine residual_and_jacobian(Tabs, Tsub, Tout, F, J, &
       G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC, &
       A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs, &
       sigma, mdot, Cp, U, A_hx, P_PV, &
       Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV)

  implicit none
  Double Precision, intent(in) :: Tabs, Tsub, Tout
  Double Precision, intent(out):: F(3), J(3,3)
  Double Precision, intent(in) :: G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC
  Double Precision, intent(in) :: A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs
  Double Precision, intent(in) :: eps_PV, eps_abs, sigma
  Double Precision, intent(in) :: mdot, Cp, U, A_hx, P_PV
  Double Precision, intent(in) :: Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV

  Double Precision :: q_solar_abs, q_solar_ptc
  Double Precision :: rad_abs, rad_pv, conv_pv, conv_abs
  Double Precision :: DT, eps_log, denom, Lmtd, Rcond_sum

  Rcond_sum = Rcond_abs + Rcond_sub
  q_solar_abs = G_GHI*A_abs*alpha_abs
  q_solar_ptc = G_DNI*A_PV*alpha_PV*CR_PTC*0.83d0*0.86d0
  rad_abs = eps_abs*sigma*A_abs*(Tabs**4 - Tsky**4)
  rad_pv  = eps_PV *sigma*A_PV*(TPV**4  - TPTC**4)
  conv_pv = (TPV - Ta)/Rconv_PV
  conv_abs= (Tabs - Ta)/Rconv_abs

  DT = Tout - Tin
  eps_log = 1.0d-9
  denom = max( (Tsub - Tout), eps_log )
  Lmtd = log( max( (Tsub - Tin), eps_log ) / denom )

  F(1) = q_solar_abs + q_solar_ptc - P_PV - rad_abs - rad_pv - conv_pv - conv_abs - mdot*Cp*DT
  F(2) = q_solar_abs - rad_abs - conv_abs - (Tabs - Tsub)/Rcond_sum
  F(3) = mdot*Cp*DT - U*A_hx * ( DT / Lmtd )

  J(1,1) = -4.0d0*eps_abs*sigma*A_abs*(Tabs**3) - 1.0d0/Rconv_abs
  J(1,2) = 0.0d0
  J(1,3) = - mdot*Cp

  J(2,1) = -4.0d0*eps_abs*sigma*A_abs*(Tabs**3) - 1.0d0/Rconv_abs - 1.0d0/Rcond_sum
  J(2,2) =  1.0d0/Rcond_sum
  J(2,3) =  0.0d0

  J(3,1) = 0.0d0
  J(3,2) = U*A_hx * DT / (Lmtd*Lmtd) * ( 1.0d0/(max(Tsub - Tin,eps_log)) - 1.0d0/(denom) )
  J(3,3) = mdot*Cp - U*A_hx * ( ( Lmtd - DT/(denom) ) / (Lmtd*Lmtd) )

end subroutine residual_and_jacobian


!===============================================================
subroutine solve_3x3(A, b, x)
  implicit none
  Double Precision, intent(in)  :: A(3,3), b(3)
  Double Precision, intent(out) :: x(3)
  Double Precision :: M(3,3), det

  M = A
  det =  M(1,1)*(M(2,2)*M(3,3)-M(2,3)*M(3,2))  &
       - M(1,2)*(M(2,1)*M(3,3)-M(2,3)*M(3,1))  &
       + M(1,3)*(M(2,1)*M(3,2)-M(2,2)*M(3,1))

  if (abs(det) < 1.0d-16) then
     print *, "Jacobiano casi singular."
  end if

  x(1) = ( b(1)*(M(2,2)*M(3,3)-M(2,3)*M(3,2))  &
         - b(2)*(M(1,2)*M(3,3)-M(1,3)*M(3,2))  &
         + b(3)*(M(1,2)*M(2,3)-M(1,3)*M(2,2)) ) / det

  x(2) = ( - b(1)*(M(2,1)*M(3,3)-M(2,3)*M(3,1))  &
           + b(2)*(M(1,1)*M(3,3)-M(1,3)*M(3,1))  &
           - b(3)*(M(1,1)*M(2,3)-M(1,3)*M(2,1)) ) / det

  x(3) = ( b(1)*(M(2,1)*M(3,2)-M(2,2)*M(3,1))  &
         - b(2)*(M(1,1)*M(3,2)-M(1,2)*M(3,1))  &
         + b(3)*(M(1,1)*M(2,2)-M(1,2)*M(2,1)) ) / det

end subroutine solve_3x3