11

In graphing a certain function using pgfplots (link here) I noticed a distortion in the curve that at first I assumed must have been a point of inflexion - but on investigating the maths, found that was not the case, the curve being very close to a straight line around that point. I have narrowed down the problem to the MWE below:

enter image description here

\documentclass[a4paper]{article}

\usepackage{tikz}
\usepackage[margin=0.3in]{geometry}
\usepackage{pgfplots}
\pgfplotsset{width=12cm,compat=1.16}

\usetikzlibrary{math}

\newcommand\alphazero{23.44} % Earth axial tilt
\newcommand\latitude{52} % latitude in degrees

\begin{document}

\begin{tikzpicture}[
    declare function={
        sunrise(\d,\x) = acos( sqrt( cos(\d)^2 - (sin(\alphazero)*cos(\x))^2 ) / cos(\d) );
    }
]

\begin{axis}[
    axis lines=left,
    align=center,
    grid=both,
    minor y tick num=4,
    title={\Large Sunrise},
    xlabel={Day of year angle $(x^{\circ})$},
    ylabel={Sunrise position $(\theta^{\circ})$},
]


% Plot (1)
% Using table calculated from Perl
\addplot [
    green
] table {perltable.dat};
\addlegendentry{(1) Table calculated from Perl};


% Plot (2) 
% Using pgfplots `\addplot expression'
\addplot expression [
    blue,
    domain=82:90,
    only marks,
    mark size=1pt,
    samples=9,
    variable=x,
]
{sunrise(\latitude, x)};
\addlegendentry{(2) Using pgfplots `\textbackslash addplot expression'};


% Plot (3)
% Plot the Tikz Math library results
\tikzmath {
    \a = sunrise(\latitude, 86);
    \b = sunrise(\latitude, 87);
    \c = sunrise(\latitude, 88);
    \d = sunrise(\latitude, 89);
    \e = sunrise(\latitude, 90);
}

\path (axis cs:82.7,0.3) node[draw,fill=white,inner sep=3pt,anchor=south west,align=left] {
    \em Tikz Math library \\
    \em results :- \\
    $\theta(86)$ = \a \\
    $\theta(87)$ = \b \\
    $\theta(88)$ = \c \\
    $\theta(89)$ = \d \\
    $\theta(90)$ = \e
};

\addplot [
    red,
    only marks,
    mark size=1pt
] coordinates {
    (86,\a) (87,\b) (88,\c) (89,\d) (90,\e)
};
\addlegendentry{(3) Tikz Math library results};

\end{axis}

\end{tikzpicture}

\end{document}


File perltable.dat:

 82.00         5.1591150521
 83.00         4.5162416308
 84.00         3.8725600710
 85.00         3.2281870745
 86.00         2.5832387643
 87.00         1.9378307813
 88.00         1.2920783808
 89.00         0.6460965291
 90.00         0.0000000000


For plot (1) I have used a Perl script to calculate the function accurately (I understand Perl computes with double-precision floating point values), and that is graphed in green. In plot (2) I use pgfplots \addplot expression, graphed in the blue dots. It can be seen there is a quite large discrepancy from the green line from x=87 onwards - this is a 50% and greater error. I also carried out another test in plot (3) using the Tikz Math library shown in the red dots, and these show a noticeable discrepancy from the blue dots - even though both plots (2) and (3) are calling the same function "sunrise".

My questions are:

  1. Is there any change to my code I can make to fix the problem of the blue dots being so far out from the green line, for example use of certain options? Or is this problem due to an inherent limitation in the accuracy of pgfplots? I am surprised the discrepancy is as large as it is. It is sufficiently big to cause us to draw a wrong mathematical conclusion about the function from the graph.

  2. Why is there such a discrepancy between the red and blue dots when they are both computed using the same function sunrise?

(In case anyone is interested this function computes the position of the sunrise as a function of latitude and day of the year).

Improve this question
0

3 Answers 3

Reset to default
6

updated

Simply to remove some completely unneeded hack of \pgfmathfloatvalueof which I had been using for the xfp and xint answers. I had been working starting with some file from initial answer now at https://tex.stackexchange.com/a/425332/4686 and only today do I read again that answer and see that later on I had removed the hack. I did not realize my starting point was only initial answer not the final version.

So now I use directly \pgfmathfloatvalueof which is not f-expandable (or was not in April 2018) but it does not need to be, pure expandability suffices for xfp or xintexpr expressions. For shortening and aesthetics I do \def\FVof(#1){\pgfmathfloatvalueof{#1}} simply to have only parentheses within the declare function expression, no braces.

You can

[x] Think a bit and use a better mathematical formula for the given range

[x] or Use xfp which achieves 16 decimal digits of precision. For this I used the technique I had developed at https://tex.stackexchange.com/a/425332/4686, but here (as I am a good boy) to plug-in xfp's \fpeval rather than \xintthefloatexpr.

[x] or Bug the xint author to implement mathematical functions, or use the available interface to program them yourselves via series in suitable ranges.

I illustrate first alternative here: (see second and third, which were added later, below)

\documentclass[a4paper]{article}

\usepackage{tikz}
\usepackage[margin=0.3in]{geometry}
\usepackage{pgfplots}
\pgfplotsset{width=12cm,compat=1.16}

\usetikzlibrary{math}

\newcommand\alphazero{23.44} % Earth axial tilt
\newcommand\latitude{52} % latitude in degrees

\begin{document}

\begin{tikzpicture}[
    % declare function={
    %     sunrise(\d,\x) = {acos( sqrt( cos(\d)^2 - (sin(\alphazero)*cos(\x))^2 ) / cos(\d) )};
    % }
    declare function={
        sunrise(\d,\x) = {asin( sin(\alphazero)*cos(\x) / cos(\d) )};
    }
]

\begin{axis}[
    axis lines=left,
    align=center,
    grid=both,
    minor y tick num=4,
    title={\Large Sunrise},
    xlabel={Day of year angle $(x^{\circ})$},
    ylabel={Sunrise position $(\theta^{\circ})$},
]


% % Plot (1)
% % Using table calculated from Perl
% \addplot [
%     green
% ] table {perltable.dat};
% \addlegendentry{(1) Table calculated from Perl};


% Plot (2) 
% Using pgfplots `\addplot expression'
\addplot expression [
    blue,
    domain=82:90,
    only marks,
    mark size=1pt,
    samples=9,
    variable=x,
]
{sunrise(\latitude, x)};
\addlegendentry{(2) Using pgfplots `\textbackslash addplot expression'};


% Plot (3)
% Plot the Tikz Math library results
\tikzmath {
    \a = sunrise(\latitude, 86);
    \b = sunrise(\latitude, 87);
    \c = sunrise(\latitude, 88);
    \d = sunrise(\latitude, 89);
    \e = sunrise(\latitude, 90);
}

\path (axis cs:82.7,0.3) node[draw,fill=white,inner sep=3pt,anchor=south west,align=left] {
    \em Tikz Math library \\
    \em results :- \\
    $\theta(86)$ = \a \\
    $\theta(87)$ = \b \\
    $\theta(88)$ = \c \\
    $\theta(89)$ = \d \\
    $\theta(90)$ = \e
};

\addplot [
    red,
    only marks,
    mark size=1pt
] coordinates {
    (86,\a) (87,\b) (88,\c) (89,\d) (90,\e)
};
\addlegendentry{(3) Tikz Math library results};

\end{axis}

\end{tikzpicture}

\end{document}

Produces:

enter image description here

For comparison I used Maple with 20 digits precision and the original (numerically not good) mathematical formula:

> sunrise(52.,86.);                                                            
                                                            2
180.00000000000000000 arccos((cos(0.28888888888888888889 Pi)

                                     2                               2 1/2
     - sin(0.13022222222222222222 Pi)  cos(0.47777777777777777778 Pi) )   /

    cos(0.28888888888888888889 Pi))/Pi

> evalf(%);
                             2.5832387643466301795

> evalf(sunrise(52.,87.));
                             1.9378307812905049847

> evalf(sunrise(52.,88.));
                             1.2920783807501711501

> evalf(sunrise(52.,89.));
                            0.64609652908098415201

One can see that with the better formula TikZ achieves about 3 or 4 digits accuracy which is enough for plotting.

Here is Maple with the better formula (and Digits set to 20):

> sunrise2:=(d,x)-> arcsin(sin(conv*23.44)*cos(conv*x)/cos(conv*d))/conv;
                                        sin(conv 23.44) cos(conv x)
                                 arcsin(---------------------------)
                                                cos(conv d)
           sunrise2 := (d, x) -> -----------------------------------
                                                conv

> evalf(sunrise2(52.,86.));                                              
                             2.5832387643466301725

> evalf(sunrise2(52.,87.));
                             1.9378307812905049928

> evalf(sunrise2(52.,88.));
                             1.2920783807501711638

> evalf(sunrise2(52.,89.));
                            0.64609652908098417757

With xfp used within declare function.

\documentclass[a4paper]{article}
\usepackage{tikz}
\usepackage[margin=0.3in]{geometry}
\usepackage{pgfplots}
\pgfplotsset{width=12cm,compat=1.16}

\usetikzlibrary{math}

\newcommand\alphazero{23.44} % Earth axial tilt
\newcommand\latitude{52} % latitude in degrees

\usepackage{xintkernel}
\usepackage{xfp}

\newcommand\FVof{}
\def\FVof(#1){\pgfmathfloatvalueof{#1}}

\begin{document}

\begin{tikzpicture}[
    declare function={
        sunrise(\d,\x) = {\fpeval{acosd( sqrt( cosd(\FVof(\d))^2 
            - (sind(\alphazero)*cosd(\FVof(\x)))^2 ) / cosd(\FVof(\d)) )}};
    }
  ]

\begin{axis}[
    axis lines=left,
    align=center,
    grid=both,
    minor y tick num=4,
    title={\Large Sunrise},
    xlabel={Day of year angle $(x^{\circ})$},
    ylabel={Sunrise position $(\theta^{\circ})$},
]


% Plot (2) 
% Using pgfplots `\addplot expression' + xfp
\addplot expression [
    blue,
    domain=82:90,
    only marks,
    mark size=1pt,
    samples=90,
    variable=x,
]
{sunrise(\latitude, x)};
\addlegendentry{Using pgfplots `\textbackslash addplot expression' and xfp and
a hack to use the latter};

% Plot (3)
% Plot the Tikz Math library results
\tikzmath {
    \a = sunrise(\latitude, 86);
    \b = sunrise(\latitude, 87);
    \c = sunrise(\latitude, 88);
    \d = sunrise(\latitude, 89);
    \e = sunrise(\latitude, 90);
}

\path (axis cs:82.7,0.3) node[draw,fill=white,inner sep=3pt,anchor=south west,align=left] {
    \em xfp results : \\
    $\theta(86)$ = \a \\
    $\theta(87)$ = \b \\
    $\theta(88)$ = \c \\
    $\theta(89)$ = \d \\
    $\theta(90)$ = \e
};

\addplot [
    red,
    only marks,
    mark size=1pt
] coordinates {
    (86,\a) (87,\b) (88,\c) (89,\d) (90,\e)
};
\addlegendentry{(3) a selection of xfp results};

\end{axis}

\end{tikzpicture}

\end{document}

enter image description here

With xint, but we needed to code sin, cos using high level interface. The most delicate is probably asind which I bothered defining only for usage with a small argument as enough here. I use for this illustration the "better" formula, because that was already quite an hurdle to implement by hand at high level the needed functions...

Compare the results with those of Maple (which was set at 20 digits of precision, whereas here we target only 16 digits) given above. Not bad...

I removed this code to make room for newer one, below.

enter image description here

I have now redone acosd/asind for xint differently, using Newton iteration rather than a series. This allows to cover the whole range.

I tested both with the "BAD" formula with acosd and the "GOOD" formula with asind. On the left xfp, on the right xint.

See code comments, for doing the plots the quantities depending only on latitude are pre-computed, but the full formula (compilation is not much longer) is also given, commented out.

\documentclass[tikz, border=12pt]{standalone}
\usepackage[T1]{fontenc}
\usepackage{tikz}
%\usepackage[margin=0.3in]{geometry}
\usepackage{pgfplots}
\pgfplotsset{width=12cm,compat=1.16}

\usetikzlibrary{math}

\newcommand\alphazero{23.44} % Earth axial tilt
\newcommand\latitude{52} % latitude in degrees

\usepackage{xfp}
\usepackage{xintexpr}

\newcommand\FVof{}
\def\FVof(#1){\pgfmathfloatvalueof{#1}}

\xintverbosetrue

% WARNING WARNING WARNING WARNING
% I HAVE TESTED ***** NOTHING ***** OF THE FOLLOWING
% EXCEPT THAT IT WORKS FOR THE CURRENT APPLICATION.
% THERE MAY BE TRIVIAL ERRORS

% \xintdeffloatvar Pi      := 3.141592653589793;

\xintdeffloatvar Piover2 := 1.570796326794897;

% \xintdeffloatvar Piover4 := 0.7853981633974483;

\xintdeffloatvar Piover180 := 0.01745329251994330;

\xintdeffloatvar invPiover180 := 57.2957795130823;

% pre-compute 1/n! for n = 1, 2, ..., 20
\xintdeffloatvar invfact\xintListWithSep{, invfact}{\xintSeq{1}{20}}
                 := seq(1/i!, i = 1..20);

% should I rather use successive divisions by (2n+1)(2n), or rather
% multiplication by their precomputed inverses, in a modified Horner scheme ?
\xintdeffloatfunc sinaux(X) := 1 - X(invfact3 - X(invfact5 - X(invfact7
                                 - X(invfact9 - X(invfact11 - X(invfact13
                                 - X(invfact15 - X * invfact17)))))));

\xintdeffloatfunc cosaux(X) := 1 - X(invfact2 - X(invfact4 - X(invfact6
                                 - X(invfact8 - X(invfact10 - X(invfact12
                                 - X(invfact14 - X(invfact16 - X * invfact18))))))));

% use this only between -pi/4 and pi/4
\xintdeffloatfunc sinAA(x) := x * sinaux(sqr(x));

% use this only between -pi/4 and pi/4
\xintdeffloatfunc cosAA(x) := cosaux(sqr(x));

% Use 1 - 2 sin(x/2)^2 formula for better cos(x) in -pi/4 < x < pi/4 ?
% And use 2 sin(x/2) cos(x/2) for sine ?

% WARNING WARNING WARNING WARNING
% I HAVE TESTED ***** NOTHING ***** OF THE FOLLOWING
% EXCEPT THAT IT WORKS FOR THE CURRENT APPLICATION.
% THERE MAY BE TRIVIAL ERRORS

\xintdeffloatfunc sind_0(x) := sinAA(x * Piover180);
\xintdeffloatfunc sind_1(x) := cosAA((x - 90) * Piover180);
\xintdeffloatfunc sind_2(x) := sind_1(x);
\xintdeffloatfunc sind_3(x) := - sinAA((x - 180) * Piover180);
\xintdeffloatfunc sind_4(x) := sind_3(x);
\xintdeffloatfunc sind_5(x) := - cosAA((x - 270) * Piover180);
\xintdeffloatfunc sind_6(x) := sind_5(x);
\xintdeffloatfunc sind_7(x) := sinAA((x - 360) * Piover180);

\xintdeffloatfunc cosd_0(x) := cosAA(x * Piover180);
\xintdeffloatfunc cosd_1(x) := -sinAA((x - 90) * Piover180);
\xintdeffloatfunc cosd_2(x) := cosd_1(x);
\xintdeffloatfunc cosd_3(x) := -cosAA((x - 180) * Piover180);
\xintdeffloatfunc cosd_4(x) := cosd_3(x);
\xintdeffloatfunc cosd_5(x) := sinAA((x - 270) * Piover180);
\xintdeffloatfunc cosd_6(x) := cosd_5(x);
\xintdeffloatfunc cosd_7(x) := cosAA((x - 360) * Piover180);

% Not trying to define it as a genuine function due to original
% dispatch, perhaps I can succeed with other syntax, of course
% goal is to evaluate only right function.

\xintNewFunction{sind_}[1]{sind_\xinttheiiexpr num(#1)//45\relax(#1)}

\xintdeffloatfunc sind(x) := sind_(x/:360);

\xintNewFunction{cosd_}[1]{cosd_\xinttheiiexpr num(#1)//45\relax(#1)}

\xintdeffloatfunc cosd(x) := cosd_(x/:360);


% WARNING WARNING WARNING WARNING
% I HAVE TESTED ***** NOTHING ***** OF THE FOLLOWING
% EXCEPT THAT IT WORKS FOR THE CURRENT APPLICATION.
% THERE MAY BE TRIVIAL ERRORS

% Compute asin(x) via Newton method
% We do the computation via a fixed point termination,
% I sense this is not very good away (should use sinaux ?)
% but will stick with that for now.
% And of course it should be written at lower level for efficiency.
\xintNewFunction{asinAA}[1]{%
   % we assume 0<= t=#1<= 0.72 and want x with sin(x) - t = 0
   % start with x=t, and iterate x <- x + (t - sin(x))/cos(x)
   % deltas will always be positive (if t non-zero)
   iter(#1;%
        subs((D<=1e-16)? % Am I sure no neverending loop due to rounding ?
                         % I should do with increased precision 
                         % and round at end
             {break(@+D)}{@+D}, 
              D=(#1 - sinAA(@))/cosAA(@)),
        i=1++ % dummy iteration index, not used but needed by iter()
   )
}

% non-negative argument only
\xintdeffloatfunc asinA(t) := 
    if(t<0.72, asinAA(t), Piover2 - asinAA(sqrt(1-sqr(t))));
\xintdeffloatfunc asindA(t) := 
    if(t<0.72, asinAA(t)*invPiover180, 90 - asinAA(sqrt(1-sqr(t)))*invPiover180);

\xintdeffloatfunc asin(t) := sgn(t) * asinA(abs(t));
\xintdeffloatfunc asind(t) := sgn(t) * asindA(abs(t));

% acos(x), acosd(x) from asin() function

\xintdeffloatfunc acos(t) := Piover2 - asin(t);
\xintdeffloatfunc acosd(t):= 90 - asind(t);


% alphazero = 23.44
% latitude ("d") = 52

% The real 2-variable xint function with no pre-computation
% "BAD" formula deliberately used to stress-test math engines
% NOTICE THAT sind(alphazero) IS COMPUTED HERE AT TIME OF THIS DEFINITION
\xintdeffloatfunc sunriseBAD(d, x) :=
    acosd( sqrt( sqr(cosd(d)) - sqr(sind(\alphazero)*cosd(x))) / cosd(d) );

% "good" formula:
\xintdeffloatfunc sunrise(d,x) := asind(sind(\alphazero) * cosd(x) / cosd(d));

% THE "CHEATING" WAY (always advisable for plotting)
% AND PRECOMPUTE SOME COEFFICIENTS
% the axis of earth is less of a variable than latitude ;-)
\xintdeffloatvar alphazero := \alphazero;%
\xintdeffloatvar sindalphazero := sind(\alphazero);% sind(23.44)
\xintdeffloatvar cosdlatitude  := cosd(\latitude);
\xintdeffloatvar cosdlatitude2 := cosdlatitude**2;

\xintdeffloatfunc sunrisecheatBAD(d, x) :=
    acosd( sqrt( cosdlatitude2 - sqr(sindalphazero*cosd(x))) / cosdlatitude );

\xintdeffloatfunc sunrisecheat(d, x) :=
    asind(sindalphazero * cosd(x) / cosdlatitude);

% WE CAN ALSO CHEAT WITH XFP

\edef\cosdlatitudesquared{\fpeval{cosd(\latitude)^2}}
\edef\cosdlatitude       {\fpeval{cosd(\latitude)}}
\edef\sindalphazero      {\fpeval{sind(\alphazero)}}
\begin{document}

\begin{tikzpicture}[
   declare function={
% in real life plotting use the cheating variants:
     sunrisexint(\d,\x) = 
        {\xintfloateval{sunrisecheatBAD(\FVof(\d),\FVof(\x))}};
     sunrisexfp(\d,\x) = % cheating variant for xfp too
        {\fpeval{
      acosd( sqrt( \cosdlatitudesquared 
                   - (\sindalphazero*cosd(\FVof(\x)))^2 
                 ) / \cosdlatitude )}};
% functions doing it again and again : 
    %  sunrisexint(\d,\x) = 
    %     {\xintfloateval{sunrise(\FVof(\d),\FVof(\x))}};
    %  sunrisexfp(\d,\x) =
    %     {\fpeval{acosd( sqrt( cosd(\FVof(\d))^2
    %                      - (sind(\alphazero)*cosd(\FVof(\x)))^2 )
    %                / cosd(\FVof(\d)))}};
    }
  ]

\begin{axis}[
    axis lines=left,
    align=center,
    grid=both,
    minor y tick num=4,
    title={\Large Sunrise using «BAD» formula},
    xlabel={Day of year angle $(x^{\circ})$},
    ylabel={Sunrise position $(\theta^{\circ})$},
]


% Plot (1) 
% Using pgfplots `\addplot expression' + xfp
\addplot expression [
    green,
    domain=-90:0,
    only marks,
    mark size=1.5pt,
    samples=31,
    variable=x,
]
{sunrisexfp(\latitude, x)};
%\addlegendentry{Using pgfplots `\textbackslash addplot expression' with the
%  «acosd» formula via xfp};

% Plot (2) 
% Using pgfplots `\addplot expression' + xint
\addplot expression [
    blue,
    domain=0:90,
    only marks,
    mark size=1.5pt,
    samples=31,
    variable=x,
]
{sunrisexint(\latitude, x)};
%\addlegendentry{Using pgfplots `\textbackslash addplot expression' with the
%  «acosd» formula via xint};

% Plot (3)
% Plot the Tikz Math library results
\tikzmath {
    \W = sunrisexint(\latitude, 0);
    \Wa = sunrisexint(\latitude, 10);
    \Wb = sunrisexint(\latitude, 20);
    \Wc = sunrisexint(\latitude, 30);
    \Wd = sunrisexint(\latitude, 40);
    \Z = sunrisexint(\latitude, 50);
    \A = sunrisexint(\latitude, 81);
    \B = sunrisexint(\latitude, 82);
    \C = sunrisexint(\latitude, 83);
    \D = sunrisexint(\latitude, 84);
    \E = sunrisexint(\latitude, 85);
    \a = sunrisexint(\latitude, 86);
    \b = sunrisexint(\latitude, 87);
    \c = sunrisexint(\latitude, 88);
    \d = sunrisexint(\latitude, 89);
    \e = sunrisexint(\latitude, 90);
%
    \Wxfp = sunrisexfp(\latitude, 0);
    \Waxfp = sunrisexfp(\latitude, 10);
    \Wbxfp = sunrisexfp(\latitude, 20);
    \Wcxfp = sunrisexfp(\latitude, 30);
    \Wdxfp = sunrisexfp(\latitude, 40);
    \Zxfp = sunrisexfp(\latitude, 50);
    \Axfp = sunrisexfp(\latitude, 81);
    \Bxfp = sunrisexfp(\latitude, 82);
    \Cxfp = sunrisexfp(\latitude, 83);
    \Dxfp = sunrisexfp(\latitude, 84);
    \Exfp = sunrisexfp(\latitude, 85);
    \axfp = sunrisexfp(\latitude, 86);
    \bxfp = sunrisexfp(\latitude, 87);
    \cxfp = sunrisexfp(\latitude, 88);
    \dxfp = sunrisexfp(\latitude, 89);
    \exfp = sunrisexfp(\latitude, 90);
}

\path (axis cs:-75,0.3) node[draw,fill=white,inner sep=3pt,anchor=south
west,align=left] {
  \begin{tabular}{rcc}
    &xfp&xint\\
    $\theta(0)$ &\Wxfp&\W \\
    $\theta(10)$ &\Waxfp&\Wa \\
    $\theta(20)$ &\Wbxfp&\Wb \\
    $\theta(30)$ &\Wcxfp&\Wc \\
    $\theta(40)$ &\Wdxfp&\Wd \\
    $\theta(50)$&\Zxfp&\Z \\
    $\theta(81)$&\Axfp&\A \\
    $\theta(82)$&\Bxfp&\B \\
    $\theta(83)$&\Cxfp&\C \\
    $\theta(84)$&\Dxfp&\D \\
    $\theta(85)$&\Exfp&\E \\
    $\theta(86)$&\axfp&\a \\
    $\theta(87)$&\bxfp&\b \\
    $\theta(88)$&\cxfp&\c \\
    $\theta(89)$&\dxfp&\d \\
    $\theta(90)$&\exfp&\e
  \end{tabular}
};

\end{axis}

\end{tikzpicture}

% Recall this:
% \xintdeffloatfunc sunrisecheat(d, x):=
%        asind(sindalphazero * cosd(x) / cosdlatitude);

\begin{tikzpicture}[
   declare function={
% in real life plotting use the cheating variants:
     sunrisexint(\d,\x) = 
        {\xintfloateval{sunrisecheat(\FVof(\d),\FVof(\x))}};
     sunrisexfp(\d,\x) = % cheating variant for xfp too
        {\fpeval{asind(\sindalphazero*cosd(\FVof(\x))/\cosdlatitude)}};
    }
  ]

\begin{axis}[
    axis lines=left,
    align=center,
    grid=both,
    minor y tick num=4,
    title={\Large Sunrise using «GOOD» formula},
    xlabel={Day of year angle $(x^{\circ})$},
    ylabel={Sunrise position $(\theta^{\circ})$},
]


% Plot (1) 
% Using pgfplots `\addplot expression' + xfp
\addplot expression [
    green,
    domain=-90:0,
    only marks,
    mark size=1.5pt,
    samples=31,
    variable=x,
]
{sunrisexfp(\latitude, x)};
%\addlegendentry{Using pgfplots `\textbackslash addplot expression' with the
%  «acosd» formula via xfp};

% Plot (2) 
% Using pgfplots `\addplot expression' + xint
\addplot expression [
    blue,
    domain=0:90,
    only marks,
    mark size=1.5pt,
    samples=31,
    variable=x,
]
{sunrisexint(\latitude, x)};
%\addlegendentry{Using pgfplots `\textbackslash addplot expression' with the
%  «acosd» formula via xint};

% Plot (3)
% Plot the Tikz Math library results
\tikzmath {
    \W = sunrisexint(\latitude, 0);
    \Wa = sunrisexint(\latitude, 10);
    \Wb = sunrisexint(\latitude, 20);
    \Wc = sunrisexint(\latitude, 30);
    \Wd = sunrisexint(\latitude, 40);
    \Z = sunrisexint(\latitude, 50);
    \A = sunrisexint(\latitude, 81);
    \B = sunrisexint(\latitude, 82);
    \C = sunrisexint(\latitude, 83);
    \D = sunrisexint(\latitude, 84);
    \E = sunrisexint(\latitude, 85);
    \a = sunrisexint(\latitude, 86);
    \b = sunrisexint(\latitude, 87);
    \c = sunrisexint(\latitude, 88);
    \d = sunrisexint(\latitude, 89);
    \e = sunrisexint(\latitude, 90);
%
    \Wxfp = sunrisexfp(\latitude, 0);
    \Waxfp = sunrisexfp(\latitude, 10);
    \Wbxfp = sunrisexfp(\latitude, 20);
    \Wcxfp = sunrisexfp(\latitude, 30);
    \Wdxfp = sunrisexfp(\latitude, 40);
    \Zxfp = sunrisexfp(\latitude, 50);
    \Axfp = sunrisexfp(\latitude, 81);
    \Bxfp = sunrisexfp(\latitude, 82);
    \Cxfp = sunrisexfp(\latitude, 83);
    \Dxfp = sunrisexfp(\latitude, 84);
    \Exfp = sunrisexfp(\latitude, 85);
    \axfp = sunrisexfp(\latitude, 86);
    \bxfp = sunrisexfp(\latitude, 87);
    \cxfp = sunrisexfp(\latitude, 88);
    \dxfp = sunrisexfp(\latitude, 89);
    \exfp = sunrisexfp(\latitude, 90);
}

\path (axis cs:-75,0.3) node[draw,fill=white,inner sep=3pt,anchor=south
west,align=left] {
  \begin{tabular}{rcc}
    &xfp&xint\\
    $\theta(0)$ &\Wxfp&\W \\
    $\theta(10)$ &\Waxfp&\Wa \\
    $\theta(20)$ &\Wbxfp&\Wb \\
    $\theta(30)$ &\Wcxfp&\Wc \\
    $\theta(40)$ &\Wdxfp&\Wd \\
    $\theta(50)$&\Zxfp&\Z \\
    $\theta(81)$&\Axfp&\A \\
    $\theta(82)$&\Bxfp&\B \\
    $\theta(83)$&\Cxfp&\C \\
    $\theta(84)$&\Dxfp&\D \\
    $\theta(85)$&\Exfp&\E \\
    $\theta(86)$&\axfp&\a \\
    $\theta(87)$&\bxfp&\b \\
    $\theta(88)$&\cxfp&\c \\
    $\theta(89)$&\dxfp&\d \\
    $\theta(90)$&\exfp&\e
  \end{tabular}
};

\end{axis}

\end{tikzpicture}

\end{document}

enter image description here

Reminder of Maple results with "bad" formula:

> evalf(sunrise(52.,86.));                                              
                             2.5832387643466301795

> evalf(sunrise(52.,87.));
                             1.9378307812905049847

> evalf(sunrise(52.,88.));
                             1.2920783807501711501

> evalf(sunrise(52.,89.));
                            0.64609652908098415201

enter image description here

Reminder of Maple results with "good" formula:

> evalf(sunrise2(52.,86.));                                              
                             2.5832387643466301725

> evalf(sunrise2(52.,87.));
                             1.9378307812905049928

> evalf(sunrise2(52.,88.));
                             1.2920783807501711638

> evalf(sunrise2(52.,89.));
                            0.64609652908098417757

Regarding xint results, in the medium range, recall that they are not using well-thought out numerical algorithms but hand-written high level routines. In particular Newton method for asind should be done with more digits precision than the target 16 precision.

Improve this answer
3
  • this answer requires xint 1.3d, because of \xintfloateval{#1}. Else use \xintthefloatexpr #1\relax as synonym valid for xint < 1.3d.
    – user4686
    Commented Jan 15, 2019 at 18:14
  • Thank you very much for these solutions you have presented. Re xfp hack - impressive that with this one piece of code, it produces 12/13 dp accuracy in \addplot compared to Maple. Re xint solution - this is achieving very high accuracy of 15 dp's compared to Maple. Looks like some pretty heroic programming! These solutions are clearly very ingenius, but they are far beyond my current level of understanding of Tex. However they are likely to be instructive for future reference......
    – Ross Ure Anderson
    Commented Jan 16, 2019 at 15:41
  • @JosephWright, in this answer I have shown how to plug-in xfp into pgfplots ; and also I have done user level implementation of sin, cos, sind, cosd, asind, acosd showing xint has potential for this. I am tempted to release some provisory package with a more complete such approximate library of math functions, because I sense implementation at core xint level might not be done anytime soon (I first need to insert some "typing" to reduce considerably the overhead in float computations, as xint was primarily a big-integer engine).
    – user4686
    Commented Jan 16, 2019 at 20:07
5

Without knowing in detail it is most likely that the "bending" is caused because of TeX's precision. There are plenty of similar questions here on TeX.SX. One way out is using gnuplot or Lua as calculation engines.

For more details please have a look at the comments in the code.

(Of course using gnuplot you need it installed, findable in the path and you need --shell-escape enabled. To use Lua as calculation engine you need to compile with LuaLaTeX.)

% used PGFPlots v1.16
\documentclass[border=5pt]{standalone}
\usepackage{pgfplots}
    \usetikzlibrary{math}
    \pgfplotsset{
        compat=1.16,
        width=12cm,
    }
    \newcommand\alphazero{23.44} % Earth axial tilt
    \newcommand\latitude{52} % latitude in degrees
\begin{document}
\begin{tikzpicture}[
    declare function={
        sunrise(\d,\x) = acos( sqrt( cos(\d)^2 - (sin(\alphazero)*cos(\x))^2 ) / cos(\d) );
        %
        % To use `Lua` as calculation engine *no* \TeX macros are allowed
        alphazero = 23.44;
        latitude = 52;
        sunrise2(\d,\x) = acos( sqrt( cos(\d)^2 - (sin(alphazero)*cos(\x))^2 ) / cos(\d) );
    },
]
    \begin{axis}[
        axis lines=left,
        align=center,
        grid=both,
        minor y tick num=4,
        title={\Large Sunrise},
        xlabel={Day of year angle $(x^{\circ})$},
        ylabel={Sunrise position $(\theta^{\circ})$},
        % (moved common options here)
        samples=9,
        domain=82:90,
        only marks,
    ]

%        % I commented this, because you didn't provide the table data
%        \addplot [
%            green
%        ] table {perltable.dat};
%            \addlegendentry{(1) Table calculated from Perl};

        \addplot+ [blue] {sunrise(\latitude, x)};
            \addlegendentry{(2) Using pgfplots `\textbackslash addplot expression'};

        % Plot (3)
        % Plot the Tikz Math library results
        \tikzmath {
            \a = sunrise(\latitude, 86);
            \b = sunrise(\latitude, 87);
            \c = sunrise(\latitude, 88);
            \d = sunrise(\latitude, 89);
            \e = sunrise(\latitude, 90);
        }

        \path (axis cs:82.7,0.3)
            node [draw,fill=white,inner sep=3pt,anchor=south west,align=left] {
            \em Tikz Math library \\
            \em results :- \\
            $\theta(86)$ = \a \\
            $\theta(87)$ = \b \\
            $\theta(88)$ = \c \\
            $\theta(89)$ = \d \\
            $\theta(90)$ = \e
        };

        \addplot+ [red] coordinates {
            (86,\a) (87,\b) (88,\c) (89,\d) (90,\e)
        };
            \addlegendentry{(3) Tikz Math library results};

        \addplot+ [brown,mark=o,mark size=2pt] gnuplot
            {acos( sqrt( cos(\latitude*pi/180)^2 - (sin(\alphazero*pi/180)*cos(x*pi/180))^2 ) / cos(\latitude*pi/180) ) *180/pi};
%        % This would be the long version using `raw gnuplot`
%        \addplot+ [brown,mark=o,mark size=2pt] gnuplot [raw gnuplot] {
%            set angles degrees;
%            set xrange [82:90];
%            set samples 9;
%            %
%            alphazero = 23.44;
%            latitude = 52;
%            sunrise(d,x) = acos( sqrt( cos(d)**2 - (sin(alphazero)*cos(x))**2 ) / cos(d) );
%            %
%            plot sunrise(latitude,x);
%        };
            \addlegendentry{(4) Using gnuplot as calculation engine for `\textbackslash addplot'};

        \addplot+ [green] {sunrise2(latitude, x)};
            \addlegendentry{(5) Using Lua as calculation engine for `\textbackslash addplot'};

    \end{axis}
\end{tikzpicture}
\end{document}

image showing the result of above code

Improve this answer
6
  • Great work! Do you use a special key combination to produce % at the start of the line instead of respecting indentation?
    – manooooh
    Commented Jan 15, 2019 at 19:18
  • 1
    @manooooh, thank you for the credits. Yes I use a special key combination in my editor to comment blocks. When you uncomment the lines starting with a % indention should be right, right? Do you think this is wrong and I should replace the spaces by %s instead of just "inserting" the %s?
    – Stefan Pinnow
    Commented Jan 15, 2019 at 19:50
  • I am not sure what you are asking, but I was not judging you by the way you comment a line (although I would respect indentation), but simply if you used some command so that the % character stays at the beginning of the line and does not accompany the indentation. I see that it is like that.
    – manooooh
    Commented Jan 15, 2019 at 20:54
  • Thank you very much, I tried both gnuplot and lualatex. '\addplot gnuplot' works very well and all the points lie along the required straight line.......
    – Ross Ure Anderson
    Commented Jan 16, 2019 at 15:42
  • For lualatex I changed the Build->Default Compiler option in TeXstudio v2.12.10 on Windows 7 from pdflatex to lualatex and then changed to function sunrise2 as you describe, then '\addplot expression' plots accurately. However I notice this does not work for the tikzmath environment - on changing it to use sunrise2 the points are just as much off as before. I am not familiar with lualatex at all, so the gnuplot solution is fine for me. Many thanks again.
    – Ross Ure Anderson
    Commented Jan 16, 2019 at 15:42
0

I believe for computing and its accuracy it is always a good choice to employ another language. I did so with FORTRAN 90 and offered the template.f90 source code at https://github.com/LiuGangKingston/FORTRAN-CSV-TIKZ.git By using template.f90, users only need to code their specific computational work in the routine "mycomputing()", as all the rest were coded already. See an example result https://github.com/LiuGangKingston/FORTRAN-CSV-TIKZ/blob/main/examples/example01.010beams/FORTRAN_CSV_TIKZ_example01.pdf the angles were calculated accurately.

Improve this answer

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .