(* ===================================== *)
(* 1) Parámetros del sistema *)
(* ===================================== *)
ClearAll["Global`*"];
m = 1;
l = 2;
R = 1;
g = 9.81;
\[Omega] = 1.74;
T = (2 \[Pi])/\[Omega];  (* periodo del forzamiento *)
(* Tiempo total y bloque de renormalización *)
Tmax = 100 T;
\[Tau] = 
  T;                 (* renormalizar cada periodo suele funcionar \
bien *)
transientSteps = 50;        (* descartar transitorio: 50 periodos *)

(* ===================================== *)
(* 2) Rejilla de condiciones iniciales *)
(* ===================================== *)
\[Phi]min = -(\[Pi]/2); \[Phi]max = \[Pi]/2;
vmin = -5; vmax = 5;
\[CapitalDelta]\[Phi] = 
  0.01;   (* AJUSTA: 0.05 o 0.1 recomendado al inicio *)
\[CapitalDelta]v = 0.01;         (* AJUSTA *)
CI = Flatten[
   Table[{\[Phi]0, 
     v0}, {\[Phi]0, \[Phi]min, \[Phi]max, \[CapitalDelta]\[Phi]}, \
{v0, vmin, vmax, \[CapitalDelta]v}  ], 1];

(* ===================================== *)
(* 3) Paso de integración (autónomo con fase psi) *)
(* ===================================== *)
stepEvolve[state_List, pert_List, \[Tau]_?NumericQ] :=
   Module[{phi0, v0, psi0, dphi0, dv0, sol, phi1, v1, psi1, dphi1, 
    dv1},
     {phi0, v0, psi0} = state;
     {dphi0, dv0} = pert;
     sol = Quiet @ Check[
           NDSolveValue[
             { \[Phi]'[t] == v[t],
               
        v'[t] == 1/
          l (R \[Omega]^2 Cos[\[Phi][t] - \[Psi][t]] - 
            g Sin[\[Phi][t]]),
               \[Psi]'[t] == \[Omega],
               (* sistema variacional (Benettin): d' = J d *)
               d\[Phi]'[t] == dv[t],
               
        dv'[t] == -(1/
           l) (R \[Omega]^2 Sin[\[Phi][t] - \[Psi][t]] + 
            g Cos[\[Phi][t]]) d\[Phi][t],
               (* CI *)
               \[Phi][0] == phi0, v[0] == v0, \[Psi][0] == psi0,
               d\[Phi][0] == dphi0, dv[0] == dv0 }, {\[Phi][\[Tau]], 
        v[\[Tau]], \[Psi][\[Tau]], d\[Phi][\[Tau]], dv[\[Tau]]}, {t, 
        0, \[Tau]},
             MaxStepSize -> 0.001, AccuracyGoal -> 12, 
       PrecisionGoal -> 12,
             Method -> {"StiffnessSwitching"}, 
       MaxSteps -> Infinity], $Failed];
     If[sol === $Failed, Return[$Failed]];
     {phi1, v1, psi1, dphi1, dv1} = sol;
     phi1 = Mod[phi1, 2 \[Pi], -\[Pi]];
     {phi1, v1, psi1, dphi1, dv1}];

(* ===================================== *)
(* 4) MLE (Benettin con renormalización periódica) *)
(* ===================================== *)
MLE[ci_List, \[Delta]0_ : 10^-8, Tmax_ : Tmax, \[Tau]_ : \[Tau],
      transientSteps_ : transientSteps] :=
   Module[{steps, state, pert, sumLog = 0., res, phi1, v1, psi1, 
    dphi1, dv1, dist, i}, steps = Floor[Tmax/\[Tau]];
     If[steps <= transientSteps, Return[Indeterminate]];
     (* estado inicial: psi(0)=0 fija la fase del forzamiento *)
     state = {ci[[1]], ci[[2]], 0.};
     (* perturbación inicial aleatoria normalizada *)
     pert = \[Delta]0 Normalize[RandomReal[{-1, 1}, 2]];
     For[i = 1, i <= steps, i++,
       res = stepEvolve[state, pert, \[Tau]];
       If[res === $Failed, Return[Indeterminate]];
       {phi1, v1, psi1, dphi1, dv1} = res;
       dist = Sqrt[dphi1^2 + dv1^2];
       If[(! NumericQ[dist]) || dist <= 0, 
     dist = \[Delta]0; {dphi1, dv1} = pert;];
       If[i > transientSteps, sumLog += Log[dist/\[Delta]0]];
       (* renormalización *)
       pert = (\[Delta]0/dist) {dphi1, dv1};
       state = {phi1, v1, psi1}; ];
     sumLog/((steps - transientSteps) \[Tau])];

(* ===================================== *)
(* 5) Cálculo paralelo del mapa *)
(* ===================================== *)
LaunchKernels[];
mleData = 
  ParallelMap[
   Function[{c}, 
    Module[{lam = MLE[c, 10^-8, Tmax, \[Tau], transientSteps]}, {c[[
       1]], c[[2]], lam}]], CI];
(* limpiar fallos *)
mleDataClean = 
  Select[mleData, 
   NumericQ[#[[3]]] && 
     FreeQ[#[[3]], Indeterminate | ComplexInfinity | Infinity] &];
{\[Lambda]min, \[Lambda]max} = MinMax[mleDataClean[[All, 3]]];

(* ===================================== *)
(* 6) Mapa de calor *)
(* ===================================== *)
ListDensityPlot[ mleDataClean, PlotRange -> Full, ImageSize -> 450, 
 AspectRatio -> 1, 
 FrameLabel -> {Style[
    "\phi", 16], 
   Style["\dot\phi", 16]}, LabelStyle -> Directive[14], ColorFunction -> "Rainbow",
  PlotLegends -> BarLegend[{"Rainbow", {\[Lambda]min, \[Lambda]max}}]]