3
$\begingroup$

I have written a Manipulate function that shows the changing photoperiod (number of hours of daylight per 24-hour day) to which a latitudinal migrant is exposed before, during, and after each of its annual migrations. The user can dial in the two latitudes and the Julian (ordinal) dates for the start and end of each of the two annual migrations. I would now like it to display the total daylight over a year to which the individual is exposed, depending on the combination of these adjustable parameters. In other words, I would like to estimate the area under the dynamically manipulated curve shown below from x = 0 to x = 365. Note that I am referring to the area under the BOLD (black and red) line, not the light gray lines. Of course, I have tried the Integrate function, but it is either WAY too slow (perhaps I am asking for something that is too computationally involved to do quickly) or maybe I am doing something wrong. Here is the Manipulate function without any attempt to calculate the area under the curve:

Manipulate[

 P = ArcSin[0.39795*Cos[0.2163108 + 2*ArcTan[0.9671396*Tan[0.0086*(x - 186)]]]];

 ASL = (Sin[DL*Pi/180] + Sin[SL*Pi/180]*Sin[P])/(Cos[SL*Pi/180]*Cos[P]);
 CASL = Clip[ASL];
 ANL = (Sin[DL*Pi/180] + Sin[NL*Pi/180]*Sin[P])/(Cos[NL*Pi/180]*Cos[P]);
 CANL = Clip[ANL];
 AMNL = (Sin[DL*Pi/180] + Sin[((((NL*x) - (NL*NMB) - (SL*x) + (SL*NMB))/(NME - NMB)) + SL)
  *Pi/180]*Sin[(ArcSin[.39795*Cos[.2163108 + 2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])])
  /(Cos[((((NL* x) - (NL*NMB) - (SL*x) + (SL*NMB))/(NME - NMB)) + SL)* Pi/180]*
  Cos[(ArcSin[.39795* Cos[.2163108 + 2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])]);
 CAMNL = Clip[AMNL]; 
 AMSL = (Sin[DL*Pi/180] + Sin[((((SL*x) - (SL*SMB) - (NL*x) + (NL*SMB))/(SME - SMB)) + NL)*
  Pi/180]*Sin[(ArcSin[.39795* Cos[.2163108 + 2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])])/
  (Cos[((((SL* x) - (SL*SMB) - (NL*x) + (NL*SMB))/(SME - SMB)) + NL)* Pi/180]*Cos[(ArcSin[
  .39795* Cos[.2163108 + 2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])]);
 CAMSL = Clip[AMSL];

 Show[

  Plot[24 - (24/Pi)*ArcCos[CASL], {x, 0, 365}, PlotRange -> {{0, 365}, {0, 24}}, Ticks -> 
   {{{79, "Mar 20"}, {172, "Jun 21"}, {265, "Sep 22"}, {355, "Dec 21"}}, {0, 2, 4, 6, 8, 
   10, 12, 14, 16, 18, 20, 22, 24}}, AxesOrigin -> {0, 12}, PlotStyle -> {Black, Thin}],
  Plot[24 - (24/Pi)*ArcCos[CANL], {x, 0, 365}, PlotRange -> {{0, 365}, {0, 24}}, PlotStyle 
   -> {Black, Thin}], 
  Plot[24 - (24/Pi)*ArcCos[CASL], {x, 0, NMB}, PlotRange -> {{0, 365}, {0, 24}}, PlotStyle 
   -> {Black}],
  Plot[24 - (24/Pi)*ArcCos[CAMNL], {x, NMB, NME}, PlotRange -> {{0, 365}, 
   {0, 24}}, PlotStyle -> {Red}],
  Plot[24 - (24/Pi)*ArcCos[CANL], {x, NME, SMB}, PlotRange -> {{0, 365}, {0, 24}}, 
   PlotStyle -> {Black}], 
  Plot[24 - (24/Pi)*ArcCos[CAMSL], {x, SMB, SME}, PlotRange -> {{0, 365}, {0, 24}}, 
   PlotStyle -> {Red}], 
  Plot[24 - (24/Pi)*ArcCos[CASL], {x, SME, 365}, PlotRange -> {{0, 365}, {0, 24}}, 
   PlotStyle -> {Black}]],

 {{DL, 0.8333, "Daylight Definition:"}, {0 -> "Sun Center at Horizon", 0.26667 -> 
   "Sun Top at Horizon", 0.8333 -> "Sun Top Apparent at Horizon", 
   6 -> "Civil Twilight Included", 12 -> "Nautical Twilight Included", 
   18 -> "Astronomical Twilight Included"}, Appearance -> "Open"},
 {{SL, -70, "Early/Late-Year Latitude"}, -90, 90, Appearance -> "Open"}, 
 {{NL, 70, "Mid-Year Latitude"}, -90, 90, Appearance -> "Open"}, 
 {{NMB, 93, "Early-Year Migration Begins"}, 1, NME, Appearance -> "Open"},
 {{NME, 136, "Early-Year Migration Ends"}, (Abs[NL - SL]/21) + NMB, 365, 
   Appearance -> "Open"}, 
 {{SMB, 229, "Late-Year Migration Begins"}, NME + 1, 355, Appearance -> "Open"}, 
 {{SME, 309, "Late-Year Migration Ends"}, (Abs[SL - NL]/21) + SMB, 364, 
   Appearance -> "Open"}]

photoperiod

Improve this question
$\endgroup$
2
  • $\begingroup$ Use Piecewise so your "complex" function appears as a single function, then use NIntegrate on that. $\endgroup$
    – george2079
    Oct 27, 2016 at 17:28
  • 3
    $\begingroup$ re: "not the light gray lines" there seems to be a lot of code in this question that is not at all relevant to the question. Work on simplifying it down to what is essential. $\endgroup$
    – george2079
    Oct 27, 2016 at 17:30

1 Answer 1

Reset to default
5
$\begingroup$

For this problem, since you are generating a plot, the fastest way to get the integral is to grab the data from the plot and use the trapezoidal rule:

 With[{SL = -70, NL = 70, NMB = 93, NME = 136, SMB = 229, SME = 309, 
  DL = .833}, 
 P = ArcSin[
   0.39795*Cos[0.2163108 + 2*ArcTan[0.9671396*Tan[0.0086*(x - 186)]]]];
 ASL = (Sin[DL*Pi/180] + Sin[SL*Pi/180]*Sin[P])/(Cos[SL*Pi/180]*
     Cos[P]);
 CASL = Clip[ASL];
 ANL = (Sin[DL*Pi/180] + Sin[NL*Pi/180]*Sin[P])/(Cos[NL*Pi/180]*
     Cos[P]);
 CANL = Clip[ANL];
 AMNL = (Sin[DL*Pi/180] + 
     Sin[((((NL*x) - (NL*NMB) - (SL*x) + (SL*NMB))/(NME - NMB)) + SL)*
        Pi/180]*Sin[(ArcSin[.39795*
          Cos[.2163108 + 
            2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])])/(Cos[((((NL*
               x) - (NL*NMB) - (SL*x) + (SL*NMB))/(NME - NMB)) + SL)*
       Pi/180]*Cos[(ArcSin[.39795*
         Cos[.2163108 + 2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])]);
 CAMNL = Clip[AMNL];
 AMSL = (Sin[DL*Pi/180] + 
     Sin[((((SL*x) - (SL*SMB) - (NL*x) + (NL*SMB))/(SME - SMB)) + NL)*
        Pi/180]*Sin[(ArcSin[.39795*
          Cos[.2163108 + 
            2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])])/(Cos[((((SL*
               x) - (SL*SMB) - (NL*x) + (NL*SMB))/(SME - SMB)) + NL)*
       Pi/180]*
     Cos[(ArcSin[.39795*
         Cos[.2163108 + 2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])]);
 CAMSL = Clip[AMSL];
 curve[x_] := 24 - (24/Pi) ArcCos@Piecewise[{
       {CASL, x <= NMB},
       {CAMNL, NMB < x <= NME},
       {CANL, NME < x <= SMB},
       {CAMSL, SMB < x <= SME},
       {CASL, SME < x <= 365}
       }];
 integral = 
  Quiet@NIntegrate[
    curve[x], {x, 0, 365}];(*comment this line out after testing*) 
 rr = Reap[Show[
    Plot[y = curve[x], {x, 0, NMB}, PlotRange -> {{0, 365}, {0, 24}}, 
     PlotStyle -> {Black}, EvaluationMonitor :> Sow[{x, y}]],
    Plot[y = curve[x], {x, NMB, NME}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Red}, 
     EvaluationMonitor :> Sow[{x, y}]],
    Plot[y = curve[x], {x, NME, SMB}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Black}, 
     EvaluationMonitor :> Sow[{x, y}]],
    Plot[y = curve[x], {x, SMB, SME}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Red}, 
     EvaluationMonitor :> Sow[{x, y}]],
    Plot[y = curve[x], {x, SME, 365}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Black}, 
     EvaluationMonitor :> Sow[{x, y}]]]];
  {rr[[1]], 
   Total[-Mean@#[[All, 2]] Subtract @@ #[[All, 1]] & /@ 
    Partition[Sort[rr[[2, 1]]], 2, 1]],
 integral}]

enter image description here

note in the end it wasn't really necessary to use the Piecewise form, you could reap/sow the data from your plot expressions the way they are.

EDIT: the original manipulate form with the integral added:

Manipulate[
 P = ArcSin[
   0.39795*Cos[0.2163108 + 2*ArcTan[0.9671396*Tan[0.0086*(x - 186)]]]];
 ASL = (Sin[DL*Pi/180] + Sin[SL*Pi/180]*Sin[P])/(Cos[SL*Pi/180]*
     Cos[P]);
 CASL = Clip[ASL];
 ANL = (Sin[DL*Pi/180] + Sin[NL*Pi/180]*Sin[P])/(Cos[NL*Pi/180]*
     Cos[P]);
 CANL = Clip[ANL];
 AMNL = (Sin[DL*Pi/180] + 
     Sin[((((NL*x) - (NL*NMB) - (SL*x) + (SL*NMB))/(NME - NMB)) + SL)*
        Pi/180]*Sin[(ArcSin[.39795*
          Cos[.2163108 + 
            2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])])/(Cos[((((NL*
               x) - (NL*NMB) - (SL*x) + (SL*NMB))/(NME - NMB)) + SL)*
       Pi/180]*Cos[(ArcSin[.39795*
         Cos[.2163108 + 2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])]);
 CAMNL = Clip[AMNL];
 AMSL = (Sin[DL*Pi/180] + 
     Sin[((((SL*x) - (SL*SMB) - (NL*x) + (NL*SMB))/(SME - SMB)) + NL)*
        Pi/180]*Sin[(ArcSin[.39795*
          Cos[.2163108 + 
            2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])])/(Cos[((((SL*
               x) - (SL*SMB) - (NL*x) + (NL*SMB))/(SME - SMB)) + NL)*
       Pi/180]*Cos[(ArcSin[.39795*
         Cos[.2163108 + 2*ArcTan[.9671396*Tan[.00860 (x - 186)]]]])]);
 CAMSL = Clip[AMSL];
 rr = Reap[
   Show[Plot[24 - (24/Pi)*ArcCos[CASL], {x, 0, 365}, 
     PlotRange -> {{0, 365}, {0, 24}}, 
     Ticks -> {{{79, "Mar 20"}, {172, "Jun 21"}, {265, 
         "Sep 22"}, {355, "Dec 21"}}, {0, 2, 4, 6, 8, 10, 12, 14, 16, 
        18, 20, 22, 24}}, AxesOrigin -> {0, 12}, 
     PlotStyle -> {Black, Thin}], 
    Plot[24 - (24/Pi)*ArcCos[CANL], {x, 0, 365}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Black, Thin}], 
    Plot[y = 24 - (24/Pi)*ArcCos[CASL], {x, 0, NMB}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Black}, 
     EvaluationMonitor :> Sow[{x, y}]],
    Plot[y = 24 - (24/Pi)*ArcCos[CAMNL], {x, NMB, NME}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Red}, 
     EvaluationMonitor :> Sow[{x, y}]],
    Plot[y = 24 - (24/Pi)*ArcCos[CANL], {x, NME, SMB}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Black}, 
     EvaluationMonitor :> Sow[{x, y}]],
    Plot[y = 24 - (24/Pi)*ArcCos[CAMSL], {x, SMB, SME}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Red}, 
     EvaluationMonitor :> Sow[{x, y}]],
    Plot[y = 24 - (24/Pi)*ArcCos[CASL], {x, SME, 365}, 
     PlotRange -> {{0, 365}, {0, 24}}, PlotStyle -> {Black}, 
     EvaluationMonitor :> Sow[{x, y}]]]]; 
 {rr[[1]], 
  StringTemplate["integral is `1`"]@
   Total[-Mean@#[[All, 2]] Subtract @@ #[[All, 1]] & /@ 
     Partition[Sort[rr[[2, 1]]], 2, 1]]}, {{DL, 0.8333, 
   "Daylight Definition:"}, {0 -> "Sun Center at Horizon", 
   0.26667 -> "Sun Top at Horizon", 
   0.8333 -> "Sun Top Apparent at Horizon", 
   6 -> "Civil Twilight Included", 12 -> "Nautical Twilight Included",
    18 -> "Astronomical Twilight Included"}, 
  Appearance -> "Open"}, {{SL, -70, "Early/Late-Year Latitude"}, -90, 
  90, Appearance -> "Open"}, {{NL, 70, "Mid-Year Latitude"}, -90, 90, 
  Appearance -> "Open"}, {{NMB, 93, "Early-Year Migration Begins"}, 1,
   NME, Appearance -> "Open"}, {{NME, 136, 
   "Early-Year Migration Ends"}, (Abs[NL - SL]/21) + NMB, 365, 
  Appearance -> "Open"}, {{SMB, 229, "Late-Year Migration Begins"}, 
  NME + 1, 355, 
  Appearance -> "Open"}, {{SME, 309, 
   "Late-Year Migration Ends"}, (Abs[SL - NL]/21) + SMB, 364, 
  Appearance -> "Open"}]
Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks much, george2079. Not being particularly handy at Mathematica, it is going to take me awhile to deploy your solution, but it looks like it solves my problem. Thanks again. $\endgroup$
    – kwsockman
    Oct 27, 2016 at 19:32
  • $\begingroup$ After several days of working with this, I am unable to produce the result I want. I don't doubt that the george2079 solution works, but, as a Mathematica novice, I am unable to implement it such that the output looks like the output I posted with the addition of the value for the area under the curve shown within the plot. In short, I would like to be able to drag the sliders in the Manipulate output and have that value for the area under the curve change in real time just like the curves themselves do. $\endgroup$
    – kwsockman
    Nov 1, 2016 at 18:06
  • $\begingroup$ see edit. You should take a good look and make sure I put EvaluationMonitor in all the right parts of the plot. $\endgroup$
    – george2079
    Nov 1, 2016 at 19:02
  • $\begingroup$ This did it, george2079. Thank you so much for the help. I am too much a novice at Mathematica to understand why what you did works; but it works. I did make some small changes for aesthetics. Specifically, I put the whole "{rr[[1]], StringTemplate . . . Sort[rr[[2, 1]]], 2, 1]]]}" portion within a "Column" command in order to remove the curly braces and comma and display the hours beneath the plot. I also rounded the calculation of hours to the nearest hour. The sliders do not produce very "smoothly/rapidly" generated output (it is somewhat jerky), but everything else is good. Thanks again. $\endgroup$
    – kwsockman
    Nov 3, 2016 at 18:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.