(******************  filename: nbIterateS.m (S := R/Rmin)  *******************)

<<DifferentialEquations`InterpolatingFunctionAnatomy`

Switch[S0, 1, targettime = targettimesec =
Input["Enter the targettime, that is, the number of seconds to integrate to
(SUGGESTED:  if Atype = Apos: 10^20  if  0 < Q0 << 1/2,   10^19  if  Q0 ~ 1/2
or greater;  if Atype = Azero: 10^19;  if Atype = Aneg: 10^18;  all can be
extended later): "]; 
              targettimeyr = targettimesec/(60*60*24*365),
           Smax, targettime = targettimesec =
Input["Enter the targettime, that is, the number of seconds to integrate to
(SUGGESTED:  10^179  if  Q0 >= 0,   10^18  if  Q0 < 0;  either can be extended
later -- trial and error might be needed: targettime should go up and down as Smax goes up or down): "];
              targettimeyr = targettimesec/(60*60*24*365)]

Print[StringForm[" "]]

If[S0 == 1,
   Print[StringForm["`1`(0) = \!\(S\_min\) = ``,     `3`'(0) = ``",
                     iS, N[N[S0,1000]]         ,     iS, 0         ]],
   Print[StringForm["`1`(0) = \!\(S\_max\) = ``,     `3`'(0) = ``",
                     iS, N[N[S0,1000]]         ,     iS, 0            ]] ]

Print[StringForm["targettime = ``  (seconds)  =  ``  (years)",
                  N[N[targettimesec,1000]],
                  EngineeringForm[N[targettimeyr,1000],6]      ]]
  

NDSolve`Iterate[state, targettime]

sol = NDSolve`ProcessSolutions[state]

         Ssol = S /. sol

      Sdotsol = Sdot /. sol

    Rsol[t_] := N[Rmin Ssol[t],1000]

term1sol[t_] := N[N[USimplify[c^2 A/3 Second^2],1000]] Ssol[t]

term2sol[t_] := (N[USimplify[2 c^2 B Second^2],1000]/
                             N[(3 Rmin^6),1000])/Ssol[t]^5

Sddotsol[t_] := term1sol[t] + term2sol[t]

 Hsqrsol[t_] := USimplify[c^2 (N[A/3,1000] - N[k/Rsol[t]^2,1000]
                                    - N[B/(3 Rsol[t]^6),1000])]

    Hsol[t_] := Sign[Sdotsol[t]] USimplify[Sqrt[Hsqrsol[t]]]

    Qsol[t_] := (Sddotsol[t]/Ssol[t])/(Sdotsol[t]/Ssol[t])^2


{lowtime, hightime} = InterpolatingFunctionDomain[First[S /. sol]][[1,1]]


Print[StringForm["lowtime = ``  (seconds),     hightime = ``  (seconds)",
                  N[lowtime]  ,     N[hightime]    ]]

Print[StringForm[" "]]

Print[StringForm[
  "`1`(targettime) = `2`(``) = ``,     `5`(targettime) = `6`(``) = ``",
   iS, iS, N[targettime], N[Ssol[targettime]],
   iR, iR, N[targettime], N[Rsol[targettime]] ]]

Print[StringForm[" "]]

Print[StringForm["`1`(targettime) = `2`(``) = `` = ``", 
                  iH, iH, N[targettime], N[Hsol[targettime]],
                          Convert[N[Hsol[targettime]],
                                  (Kilo Meter/Second)/(Mega Parsec)] ]]


Print[StringForm[" "]]

newtimeused = TimeUsed[] - totaltimeused

totaltimeused = totaltimeused + newtimeused

Print[StringForm["Time used this run  = `` = ``",
                  newtimeused seconds, (newtimeused/60) minutes]]

Print[StringForm["Total time used = `` = ``",
                  totaltimeused seconds, (totaltimeused/60) minutes]]

Print[StringForm["Memory in use = ``",
                  N[MemoryInUse[]/2^20,7] megabytes]]

(*************************  End of file nbIterateS.m *************************)

