MMS Book Software

This document provides data and programming examples discussed in:

Velten, Kahlen, Schmidt: Mathematical Modeling and Simulation, 2nd ed., Wiley-VCH, 2024

Note: If your browser does not display some of the interactive plots, use the Firefox browser.

1 Chapter: Principles of Mathematical Modeling

2 Chapter: Phenomenological Models

3 Chapter: Mechanistic Models I: ODE’s

4 Chapter: Mechanistic Models II: PDE’s

5 Chapter: Virtualization and Mixed Topics

6 Crashcourse Maxima Programming Examples

6.1 Maxima Quarto-Template e6a95681

6.2 Code 1aed3c8d

```{r}
#| message: false
#| class.source: numberLines
library(rim)
maxima.options(
  engine.format = "latex",
  engine.label = FALSE,
  inline.format = "latex", 
  inline.label = FALSE
  )
```

6.3 Code 86ffc53e

set_plot_option ([gnuplot_preamble, "set linetype 1 lw 6;set linetype 2 lw 6;set linetype 3 lw 6;set term pngcairo font 'Arial,24'"])$
set_draw_defaults(
   dimensions=[1200,1200]
   ,color=black
   ,line_width = 1
   ,nticks=200
   ,font= "Arial"
   ,font_size=30
   ,head_length=0.3
   ,head_angle=20
   ,line_width=2
   ,line_width=6
   )$

6.4 Code 88da5324

 erg:solve(K*(1+i)^4+G*(1+i)^3=P,i)$
 erg[1],K:1,G:1,P:3,numer;
 erg[2],K:1,G:1,P:3,numer;
 erg[3],K:1,G:1,P:3,numer;
 erg[4],K:1,G:1,P:3,numer;
 erg[2];
 s:string(erg[2]);
 s:ssubst("","i",s)$ 
 s:ssubst("SQRT","sqrt",s)$ 
 s:ssubst("KK","K",s)$
 s:ssubst("GG","G",s)$ 
 s:ssubst("PP","P",s);

Lines 1-6:

\[i=-2.658098067372152\] \[i=0.1222909134350642\] \[i=-1.248311578607314\,\left(-1\right)^{0.5}-1.232096423031456\] \[i=1.248311578607314\,\left(-1\right)^{0.5}-1.232096423031456\] \[i=\frac{\sqrt{-\frac{\sqrt{3}\,K\,\left(-\frac{\left(G+4\,K\right)^3}{K^3}+\frac{4\,\left(G+4\,K\right)\,\left(3\,G+6\,K\right)}{K\,K}-\frac{8\,\left(3\,G+4\,K\right)}{K}\right)}{2\,\sqrt{\frac{12\,K^2\,\left(\frac{P\,\sqrt{256\,K^3\,P+27\,G^4}}{2\,3^{\frac{3}{2}}\,K^3}+\frac{\frac{3\,\left(G^3+\left(6\,K-P\right)\,G^2+\left(12\,K^2+4\,P\,K\right)\,G+8\,K^3+8\,P\,K^2\right)}{K^3}-\frac{\left(3\,G+6\,K\right)\,\left(3\,G^2+12\,K\,G+12\,K^2+4\,P\,K\right)}{K\,K^2}}{6}+\frac{\left(3\,G+6\,K\right)^3}{27\,K^3}\right)^{\frac{2}{3}}+3\,G^2\,\left(\frac{P\,\sqrt{256\,K^3\,P+27\,G^4}}{2\,3^{\frac{3}{2}}\,K^3}+\frac{\frac{3\,\left(G^3+\left(6\,K-P\right)\,G^2+\left(12\,K^2+4\,P\,K\right)\,G+8\,K^3+8\,P\,K^2\right)}{K^3}-\frac{\left(3\,G+6\,K\right)\,\left(3\,G^2+12\,K\,G+12\,K^2+4\,P\,K\right)}{K\,K^2}}{6}+\frac{\left(3\,G+6\,K\right)^3}{27\,K^3}\right)^{\frac{1}{3}}-16\,K\,P}{\left(\frac{P\,\sqrt{256\,K^3\,P+27\,G^4}}{2\,3^{\frac{3}{2}}\,K^3}+\frac{\frac{3\,\left(G^3+\left(6\,K-P\right)\,G^2+\left(12\,K^2+4\,P\,K\right)\,G+8\,K^3+8\,P\,K^2\right)}{K^3}-\frac{\left(3\,G+6\,K\right)\,\left(3\,G^2+12\,K\,G+12\,K^2+4\,P\,K\right)}{K\,K^2}}{6}+\frac{\left(3\,G+6\,K\right)^3}{27\,K^3}\right)^{\frac{1}{3}}}}}-\left(\frac{P\,\sqrt{256\,K^3\,P+27\,G^4}}{2\,3^{\frac{3}{2}}\,K^3}+\frac{\frac{3\,\left(G^3+\left(6\,K-P\right)\,G^2+\left(12\,K^2+4\,P\,K\right)\,G+8\,K^3+8\,P\,K^2\right)}{K^3}-\frac{\left(3\,G+6\,K\right)\,\left(3\,G^2+12\,K\,G+12\,K^2+4\,P\,K\right)}{K\,K^2}}{6}+\frac{\left(3\,G+6\,K\right)^3}{27\,K^3}\right)^{\frac{1}{3}}+\frac{\frac{3\,G^2+12\,K\,G+12\,K^2+4\,P\,K}{3\,K^2}+\frac{\left(-1\right)\,\left(3\,G+6\,K\right)^2}{9\,K^2}}{\left(\frac{P\,\sqrt{256\,K^3\,P+27\,G^4}}{2\,3^{\frac{3}{2}}\,K^3}+\frac{\frac{3\,\left(G^3+\left(6\,K-P\right)\,G^2+\left(12\,K^2+4\,P\,K\right)\,G+8\,K^3+8\,P\,K^2\right)}{K^3}-\frac{\left(3\,G+6\,K\right)\,\left(3\,G^2+12\,K\,G+12\,K^2+4\,P\,K\right)}{K\,K^2}}{6}+\frac{\left(3\,G+6\,K\right)^3}{27\,K^3}\right)^{\frac{1}{3}}}+\frac{\left(G+4\,K\right)^2}{2\,K^2}-\frac{4\,\left(3\,G+6\,K\right)}{3\,K}}}{2}-\frac{\sqrt{\frac{12\,K^2\,\left(\frac{P\,\sqrt{256\,K^3\,P+27\,G^4}}{2\,3^{\frac{3}{2}}\,K^3}+\frac{\frac{3\,\left(G^3+\left(6\,K-P\right)\,G^2+\left(12\,K^2+4\,P\,K\right)\,G+8\,K^3+8\,P\,K^2\right)}{K^3}-\frac{\left(3\,G+6\,K\right)\,\left(3\,G^2+12\,K\,G+12\,K^2+4\,P\,K\right)}{K\,K^2}}{6}+\frac{\left(3\,G+6\,K\right)^3}{27\,K^3}\right)^{\frac{2}{3}}+3\,G^2\,\left(\frac{P\,\sqrt{256\,K^3\,P+27\,G^4}}{2\,3^{\frac{3}{2}}\,K^3}+\frac{\frac{3\,\left(G^3+\left(6\,K-P\right)\,G^2+\left(12\,K^2+4\,P\,K\right)\,G+8\,K^3+8\,P\,K^2\right)}{K^3}-\frac{\left(3\,G+6\,K\right)\,\left(3\,G^2+12\,K\,G+12\,K^2+4\,P\,K\right)}{K\,K^2}}{6}+\frac{\left(3\,G+6\,K\right)^3}{27\,K^3}\right)^{\frac{1}{3}}-16\,K\,P}{\left(\frac{P\,\sqrt{256\,K^3\,P+27\,G^4}}{2\,3^{\frac{3}{2}}\,K^3}+\frac{\frac{3\,\left(G^3+\left(6\,K-P\right)\,G^2+\left(12\,K^2+4\,P\,K\right)\,G+8\,K^3+8\,P\,K^2\right)}{K^3}-\frac{\left(3\,G+6\,K\right)\,\left(3\,G^2+12\,K\,G+12\,K^2+4\,P\,K\right)}{K\,K^2}}{6}+\frac{\left(3\,G+6\,K\right)^3}{27\,K^3}\right)^{\frac{1}{3}}}}}{4\,\sqrt{3}\,K}+\frac{\left(-1\right)\,\left(G+4\,K\right)}{4\,K}\]

Quarto setting for lines 7 and 12:

```{r}
maxima.options(engine.format = "ascii")
```

Line 7:

i = sqrt((-(sqrt(3)*K*((-(G+4*K)^3/K^3)+(4*(G+4*K)*(3*G+6*K))/(K*K)-(8*(3*G+4*\
K))/K))/(2*sqrt((12*K^2*((P*sqrt(256*K^3*P+27*G^4))/(2*3^(3/2)*K^3)+((3*(G^3+(\
6*K-P)*G^2+(12*K^2+4*P*K)*G+8*K^3+8*P*K^2))/K^3-((3*G+6*K)*(3*G^2+12*K*G+12*K^\
2+4*P*K))/(K*K^2))/6+(3*G+6*K)^3/(27*K^3))^(2/3)+3*G^2*((P*sqrt(256*K^3*P+27*G\
^4))/(2*3^(3/2)*K^3)+((3*(G^3+(6*K-P)*G^2+(12*K^2+4*P*K)*G+8*K^3+8*P*K^2))/K^3\
-((3*G+6*K)*(3*G^2+12*K*G+12*K^2+4*P*K))/(K*K^2))/6+(3*G+6*K)^3/(27*K^3))^(1/3\
)-16*K*P)/((P*sqrt(256*K^3*P+27*G^4))/(2*3^(3/2)*K^3)+((3*(G^3+(6*K-P)*G^2+(12\
*K^2+4*P*K)*G+8*K^3+8*P*K^2))/K^3-((3*G+6*K)*(3*G^2+12*K*G+12*K^2+4*P*K))/(K*K\
^2))/6+(3*G+6*K)^3/(27*K^3))^(1/3))))-((P*sqrt(256*K^3*P+27*G^4))/(2*3^(3/2)*K\
^3)+((3*(G^3+(6*K-P)*G^2+(12*K^2+4*P*K)*G+8*K^3+8*P*K^2))/K^3-((3*G+6*K)*(3*G^\
2+12*K*G+12*K^2+4*P*K))/(K*K^2))/6+(3*G+6*K)^3/(27*K^3))^(1/3)+((3*G^2+12*K*G+\
12*K^2+4*P*K)/(3*K^2)+((-1)*(3*G+6*K)^2)/(9*K^2))/((P*sqrt(256*K^3*P+27*G^4))/\
(2*3^(3/2)*K^3)+((3*(G^3+(6*K-P)*G^2+(12*K^2+4*P*K)*G+8*K^3+8*P*K^2))/K^3-((3*\
G+6*K)*(3*G^2+12*K*G+12*K^2+4*P*K))/(K*K^2))/6+(3*G+6*K)^3/(27*K^3))^(1/3)+(G+\
4*K)^2/(2*K^2)-(4*(3*G+6*K))/(3*K))/2-sqrt((12*K^2*((P*sqrt(256*K^3*P+27*G^4))\
/(2*3^(3/2)*K^3)+((3*(G^3+(6*K-P)*G^2+(12*K^2+4*P*K)*G+8*K^3+8*P*K^2))/K^3-((3\
*G+6*K)*(3*G^2+12*K*G+12*K^2+4*P*K))/(K*K^2))/6+(3*G+6*K)^3/(27*K^3))^(2/3)+3*\
G^2*((P*sqrt(256*K^3*P+27*G^4))/(2*3^(3/2)*K^3)+((3*(G^3+(6*K-P)*G^2+(12*K^2+4\
*P*K)*G+8*K^3+8*P*K^2))/K^3-((3*G+6*K)*(3*G^2+12*K*G+12*K^2+4*P*K))/(K*K^2))/6\
+(3*G+6*K)^3/(27*K^3))^(1/3)-16*K*P)/((P*sqrt(256*K^3*P+27*G^4))/(2*3^(3/2)*K^\
3)+((3*(G^3+(6*K-P)*G^2+(12*K^2+4*P*K)*G+8*K^3+8*P*K^2))/K^3-((3*G+6*K)*(3*G^2\
+12*K*G+12*K^2+4*P*K))/(K*K^2))/6+(3*G+6*K)^3/(27*K^3))^(1/3))/(4*sqrt(3)*K)+(\
(-1)*(G+4*K))/(4*K)

Line 12

 = SQRT((-(SQRT(3)*KK*((-(GG+4*KK)^3/KK^3)+(4*(GG+4*KK)*(3*GG+6*KK))/(KK*KK)-(\
8*(3*GG+4*KK))/KK))/(2*SQRT((12*KK^2*((PP*SQRT(256*KK^3*PP+27*GG^4))/(2*3^(3/2\
)*KK^3)+((3*(GG^3+(6*KK-PP)*GG^2+(12*KK^2+4*PP*KK)*GG+8*KK^3+8*PP*KK^2))/KK^3-\
((3*GG+6*KK)*(3*GG^2+12*KK*GG+12*KK^2+4*PP*KK))/(KK*KK^2))/6+(3*GG+6*KK)^3/(27\
*KK^3))^(2/3)+3*GG^2*((PP*SQRT(256*KK^3*PP+27*GG^4))/(2*3^(3/2)*KK^3)+((3*(GG^\
3+(6*KK-PP)*GG^2+(12*KK^2+4*PP*KK)*GG+8*KK^3+8*PP*KK^2))/KK^3-((3*GG+6*KK)*(3*\
GG^2+12*KK*GG+12*KK^2+4*PP*KK))/(KK*KK^2))/6+(3*GG+6*KK)^3/(27*KK^3))^(1/3)-16\
*KK*PP)/((PP*SQRT(256*KK^3*PP+27*GG^4))/(2*3^(3/2)*KK^3)+((3*(GG^3+(6*KK-PP)*G\
G^2+(12*KK^2+4*PP*KK)*GG+8*KK^3+8*PP*KK^2))/KK^3-((3*GG+6*KK)*(3*GG^2+12*KK*GG\
+12*KK^2+4*PP*KK))/(KK*KK^2))/6+(3*GG+6*KK)^3/(27*KK^3))^(1/3))))-((PP*SQRT(25\
6*KK^3*PP+27*GG^4))/(2*3^(3/2)*KK^3)+((3*(GG^3+(6*KK-PP)*GG^2+(12*KK^2+4*PP*KK\
)*GG+8*KK^3+8*PP*KK^2))/KK^3-((3*GG+6*KK)*(3*GG^2+12*KK*GG+12*KK^2+4*PP*KK))/(\
KK*KK^2))/6+(3*GG+6*KK)^3/(27*KK^3))^(1/3)+((3*GG^2+12*KK*GG+12*KK^2+4*PP*KK)/\
(3*KK^2)+((-1)*(3*GG+6*KK)^2)/(9*KK^2))/((PP*SQRT(256*KK^3*PP+27*GG^4))/(2*3^(\
3/2)*KK^3)+((3*(GG^3+(6*KK-PP)*GG^2+(12*KK^2+4*PP*KK)*GG+8*KK^3+8*PP*KK^2))/KK\
^3-((3*GG+6*KK)*(3*GG^2+12*KK*GG+12*KK^2+4*PP*KK))/(KK*KK^2))/6+(3*GG+6*KK)^3/\
(27*KK^3))^(1/3)+(GG+4*KK)^2/(2*KK^2)-(4*(3*GG+6*KK))/(3*KK))/2-SQRT((12*KK^2*\
((PP*SQRT(256*KK^3*PP+27*GG^4))/(2*3^(3/2)*KK^3)+((3*(GG^3+(6*KK-PP)*GG^2+(12*\
KK^2+4*PP*KK)*GG+8*KK^3+8*PP*KK^2))/KK^3-((3*GG+6*KK)*(3*GG^2+12*KK*GG+12*KK^2\
+4*PP*KK))/(KK*KK^2))/6+(3*GG+6*KK)^3/(27*KK^3))^(2/3)+3*GG^2*((PP*SQRT(256*KK\
^3*PP+27*GG^4))/(2*3^(3/2)*KK^3)+((3*(GG^3+(6*KK-PP)*GG^2+(12*KK^2+4*PP*KK)*GG\
+8*KK^3+8*PP*KK^2))/KK^3-((3*GG+6*KK)*(3*GG^2+12*KK*GG+12*KK^2+4*PP*KK))/(KK*K\
K^2))/6+(3*GG+6*KK)^3/(27*KK^3))^(1/3)-16*KK*PP)/((PP*SQRT(256*KK^3*PP+27*GG^4\
))/(2*3^(3/2)*KK^3)+((3*(GG^3+(6*KK-PP)*GG^2+(12*KK^2+4*PP*KK)*GG+8*KK^3+8*PP*\
KK^2))/KK^3-((3*GG+6*KK)*(3*GG^2+12*KK*GG+12*KK^2+4*PP*KK))/(KK*KK^2))/6+(3*GG\
+6*KK)^3/(27*KK^3))^(1/3))/(4*SQRT(3)*KK)+((-1)*(GG+4*KK))/(4*KK)

Quarto setting (back to default):

```{r}
maxima.options(engine.format = "latex")
```

6.5 Code ab72c45f

f(x):=cos(x)-x$
define(f1(x),diff(f(x),x))$
xold:1$
xapp:[]$
for n:1 thru 4 do (
    xnew:xold-f(xold)/f1(xold),
    xapp:append(xapp,[ev(xnew,numer)]),
    xold:xnew
)$
xapp;
find_root(cos(x)=x,x,0,1);

\[\left[ 0.7503638678402439 , 0.7391128909113617 , 0.739085133385284 , 0.7390851332151607 \right] \] \[0.7390851332151607\]

6.6 Code 1ff10ce0

f(x):=(12*x^10-360*x^9+4515*x^8-30600*x^7+120680*x^6-276192*x^5+339120*x^4-172800*x^3)/120$
define(f1(x),diff(f(x),x));
define(f2(x),diff(f1(x),x));
krit:solve(f1(x)=0);
max:[]$
min:[]$
chk:[]$
for kritp in krit do (
    beda:ev(f2(x)>0,kritp,pred), 
    bedb:ev(f2(x)<0,kritp,pred), 
     if (beda) then min:append(min,[kritp])
       elseif (bedb) then max: append(max,[kritp])
       else chk: append(chk,[kritp])
   )$
max;
min;
chk;
draw2d(explicit(f(x),x,-0.25,6.25))$

Lines 1-4:

\[f\left(x\right):=\frac{12\,x^{10}-360\,x^9+4515\,x^8-30600\,x^7+120680\,x^6-276192\,x^5+339120\,x^4-172800\,x^3}{120}\] \[f_{1}\left(x\right):=\frac{120\,x^9-3240\,x^8+36120\,x^7-214200\,x^6+724080\,x^5-1380960\,x^4+1356480\,x^3-518400\,x^2}{120}\] \[f_{2}\left(x\right):=\frac{1080\,x^8-25920\,x^7+252840\,x^6-1285200\,x^5+3620400\,x^4-5523840\,x^3+4069440\,x^2-1036800\,x}{120}\] \[\left[ x=3 , x=4 , x=2 , x=5 , x=1 , x=6 , x=0 \right] \]

Lines 16-18:

\[\left[ x=4 , x=2 \right] \] \[\left[ x=3 , x=5 , x=1 \right] \] \[\left[ x=6 , x=0 \right] \]

Line 19:

6.7 Code f409a2bd

f1(x):=x^2*(x-1)*(x-2)*(x-3)*(x-4)*(x-5)*(x-6)^2;
define(f(x),integrate(f1(x),x));

\[f_{1}\left(x\right):=x^2\,\left(x-1\right)\,\left(x-2\right)\,\left(x-3\right)\,\left(x-4\right)\,\left(x-5\right)\,\left(x-6\right)^2\] \[f\left(x\right):=\frac{12\,x^{10}-360\,x^9+4515\,x^8-30600\,x^7+120680\,x^6-276192\,x^5+339120\,x^4-172800\,x^3}{120}\]

6.8 Code c21689ff

f(x,y):=x*exp(-x^2-y^2);
fx:diff(f(x,y),x); 
fy:diff(f(x,y),y);
fxx:define(fxx(x,y),diff(diff(f(x,y),x),x)); 
fyy:define(fyy(x,y),diff(diff(f(x,y),y),y)); 
fxy:define(fxy(x,y),diff(diff(f(x,y),x),y)); 
krit:solve([diff(f(x,y),x)=0,diff(f(x,y),y)=0]);
max:[]$
min:[]$
chk:[]$
for kritp in krit do (
    beda:ev(fxx(x,y)*fyy(x,y)>fxy(x,y)^2,kritp,pred), 
    bedb:ev(fxx(x,y)>0,kritp,pred), 
    bedc:ev(fxx(x,y)<0,kritp,pred), 
     if (beda and bedb) then min:append(min,[kritp])
       elseif (beda and bedc) then max: append(max,[kritp])
       else chk: append(chk,[kritp])
   )$
max;
min;
chk;
plot3d(f(x,y),[x,-1,1],[y,-1,1])$
draw3d(explicit(f(x,y),x,-1,1,y,-1,1))$

Lines 1-7:

\[f\left(x , y\right):=x\,\exp \left(-x^2-y^2\right)\] \[e^{-y^2-x^2}-2\,x^2\,e^{-y^2-x^2}\] \[-2\,x\,y\,e^{-y^2-x^2}\] \[\textit{fxx}\left(x , y\right):=4\,x^3\,e^{-y^2-x^2}-6\,x\,e^{-y^2-x^2}\] \[\textit{fyy}\left(x , y\right):=4\,x\,y^2\,e^{-y^2-x^2}-2\,x\,e^{-y^2-x^2}\] \[\textit{fxy}\left(x , y\right):=4\,x^2\,y\,e^{-y^2-x^2}-2\,y\,e^{-y^2-x^2}\] \[\left[ \left[ y=0 , x=-\frac{1}{\sqrt{2}} \right] , \left[ y=0 , x=\frac{1}{\sqrt{2}} \right] \right] \]

Lines 19-21:

\[\left[ \left[ y=0 , x=\frac{1}{\sqrt{2}} \right] \right] \] \[\left[ \left[ y=0 , x=-\frac{1}{\sqrt{2}} \right] \right] \] \[\left[ \right] \]

Line 22-23:

6.9 Code e5007dca

f(x,y):=cos(x)*sin(y);
plot3d(f(x,y),[x,0,10],[y,0,20],[z,-1,1])$
define(fx(x,y),diff(f(x,y),x));
define(fy(x,y),diff(f(x,y),y));
erg:quad_qags(
      'quad_qags(sqrt(1+fx(x,y)^2+fy(x,y)^2),x,0,10)[1],y,0,20)[1];

Line 1:

\[f\left(x , y\right):=\cos x\,\sin y\]

Line 2:

Lines 3-6:

\[\textit{fx}\left(x , y\right):=-\sin x\,\sin y\] \[\textit{fy}\left(x , y\right):=\cos x\,\cos y\] \[244.1109878993051\]

6.10 Code cb74ba50

v:[1,2]$
w:[2,1]$
draw2d(
       vector([0,0],v)
       ,vector([0,0],w)
       );

6.11 Code 10cc72d5

v:[1,2,3]$
w:[3,2,1]$
draw3d(
       color=green,vector([0,0,0],v)
       ,color=red ,vector([0,0,0],w)
       ,color=blue,vector(v,w-v)
       );

6.12 Code e6a95681

r:1$
kreis:parametric(r*cos(phi),r*sin(phi),phi,0,2*%pi)$
v0:vector([0,0],[r*cos(phi),r*sin(phi)]),phi=0$
v45:vector([0,0],[r*cos(phi),r*sin(phi)]),phi=45*%pi/180$
v90:vector([0,0],[r*cos(phi),r*sin(phi)]),phi=90*%pi/180$
draw2d(kreis
       ,user_preamble="set size ratio -1"
       ,color=red,v0
       ,color=blue,v45
       ,color=green,v90
       );

6.13 Code cc5febc6

scene:[]$
for phi:0 thru 2*%pi step 2*%pi/12  do
  (angle:float(phi*180/%pi),
   scene: append(scene,[
   gr2d(
        title=concat("Angle: ",angle, " degrees")
        ,color=red,vector([0,0],[cos(phi),sin(phi)])
        ,color=blue, vector([0,0],[1,0])
        ,parametric(cos(x),sin(x),x,0,2*%pi)
        ,xrange=[-1,1]
        ,yrange=[-1,1]
    )]))$
  
draw(
     delay=100
     ,scene
     ,terminal=animated_gif,file_name = "cc5febc6"
    )$

```{r,eval=show1}
#| out.width: 100%
knitr::include_graphics("./cc5febc6.gif")
```

6.14 Code ac8413f0

r:1$
kugel:parametric_surface(
 r*sin(theta)*cos(phi)
,r*sin(theta)*sin(phi)
,r*cos(theta)
,theta,0,%pi
,phi,0,2*%pi
 )$
draw3d(
 proportional_axes=xyz
,enhanced3d= true
,palette = [blue, green,yellow,orange,red]
,wired_surface = true
,kugel
);

6.15 Code ec6cea18

scene:[]$
for r:0 thru 1.7 step 0.1  do
  (x:printf(false, "~,1f",round(r*10)/10),
   kugel:parametric_surface(
    r*sin(theta)*cos(phi)
   ,r*sin(theta)*sin(phi)
   ,r*cos(theta)
   ,theta,0,%pi
   ,phi,0,2*%pi
   ),
   scene: append(scene,[
   gr3d(
        title=concat("r: ",x)
        ,proportional_axes=xyz
        ,xrange = [-1,1]
        ,yrange = [-1,1]
        ,zrange = [-1,1]
        ,enhanced3d= true
        ,palette = [blue, green,yellow,orange,red]
        ,wired_surface = true
        ,axis_3d = false
        ,colorbox=false
        ,kugel
    )]))$
  
draw(
     delay=100
     ,scene
     ,terminal=animated_gif,file_name = "cc5febcx"
    )$

```{r,eval=show1}
#| out.width: 100%
knitr::include_graphics("./cc5febcx.gif")
```

6.16 Code 1fbdd006

r:1$
scene:[]$
for phiv:0 thru 2*%pi step 2*%pi/10  do
  (x:printf(false, "~,1f",round(phiv*10)/10),
   kugel:parametric_surface(
    r*sin(theta)*cos(phi)
   ,r*sin(theta)*sin(phi)
   ,r*cos(theta)
   ,theta,0,%pi
   ,phi,0,phiv
   ),
   scene: append(scene,[
   gr3d(
        title=concat("phi: ",x)
        ,proportional_axes=xyz
        ,xrange = [-1,1]
        ,yrange = [-1,1]
        ,zrange = [-1,1]
        ,enhanced3d= true
        ,palette = [blue, green,yellow,orange,red]
        ,wired_surface = true
        ,axis_3d = false
        ,colorbox=false
        ,kugel
    )]))$
  
draw(
     delay=100
     ,scene
     ,terminal=animated_gif,file_name = "1fbdd006"
    )$

```{r,eval=show1}
#| out.width: 100%
knitr::include_graphics("./1fbdd006.gif")
```

6.17 Code 27ce7837

ode1:'diff(T(t),t)=r*(T[b]-T(t));
desolve(ode1,T(t));

\[\frac{d}{d\,t}\,T\left(t\right)=r\,\left(T_{b}-T\left(t\right)\right)\] \[T\left(t\right)=\left(T\left(0\right)-T_{b}\right)\,e^ {- r\,t }+T_{b}\]

6.18 Code 8db96755

ode1:'diff(T(t),t)=r*(T[b]-T(t));
atvalue(T(t),t=0,T0)$
desolve(ode1,T(t));

\[\frac{d}{d\,t}\,T\left(t\right)=r\,\left(T_{b}-T\left(t\right)\right)\] \[T\left(t\right)=\left(T_{0}-T_{b}\right)\,e^ {- r\,t }+T_{b}\]

6.19 Code b06af1b6

params:[r=0.281,T0=32.2,T[b]=37.2]$
ode1:'diff(T(t),t)=r*(T[b]-T(t))$
atvalue(T(t),t=0,T0)$
desolve(ode1,T(t))$
define(T(t),rhs(%)),params;
plot2d(T(t),[t,0,18],[xlabel,"t [min]"],[ylabel,"Temperatur [Grad]"])$

\[T\left(t\right):=37.2-5.0\,e^ {- 0.281\,t }\]

6.20 Code 92c87653

fever.csv

data:read_nested_list("./fever.csv")$
data:rest(data,1)$
data;
plot2d([T(t),[discrete,data]]
       ,[t,0,18]
       ,[style,[lines],[points]]
       ,[xlabel,"t [min]"]
       ,[ylabel,"Temperatur [Grad]"]
       ,[legend, false]
       )$

Lines 1-3:

\[\left[ \left[ 0 , 32.2 \right] , \left[ 2 , 34.5 \right] , \left[ 4 , 35.5 \right] , \left[ 6 , 36.3 \right] , \left[ 8 , 36.6 \right] , \left[ 10 , 36.9 \right] , \left[ 12 , 37 \right] , \left[ 14 , 37.1 \right] , \left[ 18 , 37.2 \right] \right] \]

Lines 4-10:

6.21 Code dc0e180b

params:[r=0.281,T0=32.2,Tb=37.2]$
sol:rk(r*(Tb-T),T,T0,[t,0,18,2]),params$
sol;
plot2d([discrete,sol]
       ,[x,0,18]
       ,[style,points]
       ,[xlabel,"t [min]"]
       ,[ylabel,"Temperatur [Grad]"]
       )$

\[\left[ \left[ 0.0 , 32.2 \right] , \left[ 2.0 , 34.34752747493 \right] , \left[ 4.0 , 35.57268009874416 \right] , \left[ 6.0 , 36.27162293843361 \right] , \left[ 8.0 , 36.67036598779533 \right] , \left[ 10.0 , 36.89784670636872 \right] , \left[ 12.0 , 37.02762320631147 \right] , \left[ 14.0 , 37.10165998640876 \right] , \left[ 16.0 , 37.14389756262319 \right] , \left[ 18.0 , 37.16799386775864 \right] \right] \]

6.22 Code 3240b9e8

room.csv

ode1:'diff(Ts(t),t)=rsi*(Ti(t)-Ts(t));
ode2:'diff(Ti(t),t)=ria*(Ta-Ti(t));
atvalue(Ts(t),t=0,Ts0)$
atvalue(Ti(t),t=0,Ti0)$
sol:desolve([ode1,ode2],[Ts(t),Ti(t)]);
define(Ts(t),ev(Ts(t),sol));
define(Ti(t),ev(Ti(t),sol));
params:[rsi=0.18,ria=0.15,Ta=21,Ts0=18.5,Ti0=17]$
define(Ts(t),ev(Ts(t),sol,params));
define(Ti(t),ev(Ti(t),sol,params));
data:read_nested_list("./b/room.csv")$
data:rest(data,1)$
set_plot_option ([gnuplot_preamble, "set key bottom right;set linetype 1 lw 6;set linetype 2 lw 6;set linetype 3 lw 6;set term pngcairo font 'Arial,24'"])$
plot2d([Ts(t),Ti(t),[discrete,data]]
       ,[t,0,60]
       ,[style,[lines],[lines],[points]]
       ,[xlabel,"t [min]"]
       ,[ylabel,"Temperatur [Grad]"]
       ,[legend,"Sensor temperature","Inside temperature","Data"]
       )$

Lines 1-5:

\[\frac{d}{d\,t}\,\textit{Ts}\left(t\right)=\textit{rsi}\,\left(\textit{Ti}\left(t\right)-\textit{Ts}\left(t\right)\right)\] \[\frac{d}{d\,t}\,\textit{Ti}\left(t\right)=\textit{ria}\,\left(\textit{Ta}-\textit{Ti}\left(t\right)\right)\] \[\left[ \textit{Ts}\left(t\right)=\frac{\left(\left(\textit{Ts}_{0}-\textit{Ti}_{0}\right)\,\textit{rsi}+\left(\textit{Ta}-\textit{Ts}_{0}\right)\,\textit{ria}\right)\,e^ {- \textit{rsi}\,t }}{\textit{rsi}-\textit{ria}}+\frac{\left(\textit{Ti}_{0}-\textit{Ta}\right)\,\textit{rsi}\,e^ {- \textit{ria}\,t }}{\textit{rsi}-\textit{ria}}+\textit{Ta} , \textit{Ti}\left(t\right)=\left(\textit{Ti}_{0}-\textit{Ta}\right)\,e^ {- \textit{ria}\,t }+\textit{Ta} \right] \]

Lines 6-7:

\[\textit{Ts}\left(t\right):=\frac{\left(\left(\textit{Ts}_{0}-\textit{Ti}_{0}\right)\,\textit{rsi}+\left(\textit{Ta}-\textit{Ts}_{0}\right)\,\textit{ria}\right)\,e^ {- \textit{rsi}\,t }}{\textit{rsi}-\textit{ria}}+\frac{\left(\textit{Ti}_{0}-\textit{Ta}\right)\,\textit{rsi}\,e^ {- \textit{ria}\,t }}{\textit{rsi}-\textit{ria}}+\textit{Ta}\] \[\textit{Ti}\left(t\right):=\left(\textit{Ti}_{0}-\textit{Ta}\right)\,e^ {- \textit{ria}\,t }+\textit{Ta}\]

Lines 8-10:

\[\textit{Ts}\left(t\right):=-24.0\,e^ {- 0.15\,t }+21.5\,e^ {- 0.18\,t }+21\] \[\textit{Ti}\left(t\right):=21-4\,e^ {- 0.15\,t }\]

Lines 11-20:

7 Crashcourse R Programming Examples

7.1 Preliminaries

7.1.1 Course data 9411def5

7.2 Basic workflow

7.2.1 Code c55d0a7b

data/weather.csv:
da,p1,p2,p3,pd,tmax,tmin,emin,t1,t2,t3,rm1,rm2,rm3,rmd,w1,w2,w3,ws1,ws2,ws3,wsd ,cc1,cc2,cc3,ccd,sd,vr1,vr2,vr3,pr1,prt1,pr2,prt2,pr3,prt3,pd,pdt,sh,wsdmax
01.01.1961,1001.4,1001.5,1001.8,1001.6,2.8,1.4,0.0,1.5,2.6,2.4,95,90,87,91,C,C,C,0,0,0,0.0,8,8,8,8.0,0,0.5 bis 1 km,2 bis 4 km,1 bis 2 km,0.0,-,0,-,0,-,0.0,-,0,3.5
02.01.1961,996.7,992.0,985.8,991.5,4.1,0.9,-1.1,3.1,2.6,3.7,77,94,95,89,S,S,NO,2,1,2,1.7,7,8,8,7.7,0,10 bis 20 km,4 bis 10 km,2 bis 4 km,0.0,-,0.7,Regen,1.7,Regen,7.9,Regen,0,6.9
03.01.1961,973.5,980.6,981.4,978.5,8.6,2.7,3.0,8.4,4.9,3.5,78,83,77,79,S,S,SO,5,1,1,2.3,8,7,7,7.3,0,4 bis 10 km,20 bis 50 km,10 bis 20 km,5.5,Regen,1.4,Regen,0,-,1.4,Regen,0,18.4
04.01.1961,980.4,982.6,986.6,983.2,5.5,1.5,-1.3,3.4,5.4,2.6,84,85,92,87,S,S,S,1,2,2,1.7,7,6,7,6.7,0.3,4 bis 10 km,20 bis 50 km,4 bis 10 km,0.0,Regen,0.6,Regen,4.3,Regen,5.0,Regen,0,14.1
05.01.1961,993.6,1001.4,1003.9,999.6,6.4,1.7,-0.4,4.4,5.9,2,77,65,82,75,SW,SW,S,3,3,2,2.7,6,4,7,5.7,2.5,4 bis 10 km,20 bis 50 km,4 bis 10 km,0.1,Regen,0,Regen,0,-,0.0,Regen,0,15
...

d=read.csv("data/weather.csv")
d[1:5,]
class(d$da)
d$da1=as.Date(d$da,"%d.%m.%Y")
class(d$da1)
d[1:5,]
save(d,file="data/c55d0a7b.RData")

d=read.csv("data/weather.csv")
d[1:5,]

          da     p1     p2     p3     pd tmax tmin emin  t1  t2  t3 rm1 rm2 rm3
1 01.01.1961 1001.4 1001.5 1001.8 1001.6  2.8  1.4  0.0 1.5 2.6 2.4  95  90  87
2 02.01.1961  996.7  992.0  985.8  991.5  4.1  0.9 -1.1 3.1 2.6 3.7  77  94  95
3 03.01.1961  973.5  980.6  981.4  978.5  8.6  2.7  3.0 8.4 4.9 3.5  78  83  77
4 04.01.1961  980.4  982.6  986.6  983.2  5.5  1.5 -1.3 3.4 5.4 2.6  84  85  92
5 05.01.1961  993.6 1001.4 1003.9  999.6  6.4  1.7 -0.4 4.4 5.9 2.0  77  65  82
  rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd  sd          vr1          vr2
1  91  C  C  C   0   0   0 0.0   8   8   8 8.0 0.0 0.5 bis 1 km   2 bis 4 km
2  89  S  S NO   2   1   2 1.7   7   8   8 7.7 0.0 10 bis 20 km  4 bis 10 km
3  79  S  S SO   5   1   1 2.3   8   7   7 7.3 0.0  4 bis 10 km 20 bis 50 km
4  87  S  S  S   1   2   2 1.7   7   6   7 6.7 0.3  4 bis 10 km 20 bis 50 km
5  75 SW SW  S   3   3   2 2.7   6   4   7 5.7 2.5  4 bis 10 km 20 bis 50 km
           vr3 pr1  prt1 pr2  prt2 pr3  prt3 pd.1   pdt sh wsdmax
1   1 bis 2 km 0.0     - 0.0     - 0.0     -  0.0     -  0    3.5
2   2 bis 4 km 0.0     - 0.7 Regen 1.7 Regen  7.9 Regen  0    6.9
3 10 bis 20 km 5.5 Regen 1.4 Regen 0.0     -  1.4 Regen  0   18.4
4  4 bis 10 km 0.0 Regen 0.6 Regen 4.3 Regen  5.0 Regen  0   14.1
5  4 bis 10 km 0.1 Regen 0.0 Regen 0.0     -  0.0 Regen  0   15.0

class(d$da)

[1] "character"

d$da1=as.Date(d$da,"%d.%m.%Y")
class(d$da1)

[1] "Date"

d[1:5,]

          da     p1     p2     p3     pd tmax tmin emin  t1  t2  t3 rm1 rm2 rm3
1 01.01.1961 1001.4 1001.5 1001.8 1001.6  2.8  1.4  0.0 1.5 2.6 2.4  95  90  87
2 02.01.1961  996.7  992.0  985.8  991.5  4.1  0.9 -1.1 3.1 2.6 3.7  77  94  95
3 03.01.1961  973.5  980.6  981.4  978.5  8.6  2.7  3.0 8.4 4.9 3.5  78  83  77
4 04.01.1961  980.4  982.6  986.6  983.2  5.5  1.5 -1.3 3.4 5.4 2.6  84  85  92
5 05.01.1961  993.6 1001.4 1003.9  999.6  6.4  1.7 -0.4 4.4 5.9 2.0  77  65  82
  rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd  sd          vr1          vr2
1  91  C  C  C   0   0   0 0.0   8   8   8 8.0 0.0 0.5 bis 1 km   2 bis 4 km
2  89  S  S NO   2   1   2 1.7   7   8   8 7.7 0.0 10 bis 20 km  4 bis 10 km
3  79  S  S SO   5   1   1 2.3   8   7   7 7.3 0.0  4 bis 10 km 20 bis 50 km
4  87  S  S  S   1   2   2 1.7   7   6   7 6.7 0.3  4 bis 10 km 20 bis 50 km
5  75 SW SW  S   3   3   2 2.7   6   4   7 5.7 2.5  4 bis 10 km 20 bis 50 km
           vr3 pr1  prt1 pr2  prt2 pr3  prt3 pd.1   pdt sh wsdmax        da1
1   1 bis 2 km 0.0     - 0.0     - 0.0     -  0.0     -  0    3.5 1961-01-01
2   2 bis 4 km 0.0     - 0.7 Regen 1.7 Regen  7.9 Regen  0    6.9 1961-01-02
3 10 bis 20 km 5.5 Regen 1.4 Regen 0.0     -  1.4 Regen  0   18.4 1961-01-03
4  4 bis 10 km 0.0 Regen 0.6 Regen 4.3 Regen  5.0 Regen  0   14.1 1961-01-04
5  4 bis 10 km 0.1 Regen 0.0 Regen 0.0     -  0.0 Regen  0   15.0 1961-01-05

save(d,file="data/c55d0a7b.RData")

7.2.2 Code 907b39d3

data/weather2.csv:
da;p1;p2;p3;pd;tmax;tmin;emin;t1;t2;t3;rm1;rm2;rm3;rmd;w1;w2;w3;ws1;ws2;ws3;wsd ;cc1;cc2;cc3;ccd;sd;vr1;vr2;vr3;pr1;prt1;pr2;prt2;pr3;prt3;pd;pdt;sh;wsdmax
01.01.1961;1001,4;1001,5;1001,8;1001,6;2,8;1,4;0,0;1,5;2,6;2,4;95;90;87;91;C;C;C;0;0;0;0,0;8;8;8;8,0;0;0.5 bis 1 km;2 bis 4 km;1 bis 2 km;0,0;-;0;-;0;-;0,0;-;0;3,5
02.01.1961;996,7;992,0;985,8;991,5;4,1;0,9;-1,1;3,1;2,6;3,7;77;94;95;89;S;S;NO;2;1;2;1,7;7;8;8;7,7;0;10 bis 20 km;4 bis 10 km;2 bis 4 km;0,0;-;0,7;Regen;1,7;Regen;7,9;Regen;0;6,9
03.01.1961;973,5;980,6;981,4;978,5;8,6;2,7;3,0;8,4;4,9;3,5;78;83;77;79;S;S;SO;5;1;1;2,3;8;7;7;7,3;0;4 bis 10 km;20 bis 50 km;10 bis 20 km;5,5;Regen;1,4;Regen;0;-;1,4;Regen;0;18,4
04.01.1961;980,4;982,6;986,6;983,2;5,5;1,5;-1,3;3,4;5,4;2,6;84;85;92;87;S;S;S;1;2;2;1,7;7;6;7;6,7;0,3;4 bis 10 km;20 bis 50 km;4 bis 10 km;0,0;Regen;0,6;Regen;4,3;Regen;5,0;Regen;0;14,1
05.01.1961;993,6;1001,4;1003,9;999,6;6,4;1,7;-0,4;4,4;5,9;2;77;65;82;75;SW;SW;S;3;3;2;2,7;6;4;7;5,7;2,5;4 bis 10 km;20 bis 50 km;4 bis 10 km;0,1;Regen;0;Regen;0;-;0,0;Regen;0;15
...

d=read.csv2("data/weather2.csv")
d[1:5,]

          da     p1     p2     p3     pd tmax tmin emin  t1  t2  t3 rm1 rm2 rm3
1 01.01.1961 1001.4 1001.5 1001.8 1001.6  2.8  1.4  0.0 1.5 2.6 2.4  95  90  87
2 02.01.1961  996.7  992.0  985.8  991.5  4.1  0.9 -1.1 3.1 2.6 3.7  77  94  95
3 03.01.1961  973.5  980.6  981.4  978.5  8.6  2.7  3.0 8.4 4.9 3.5  78  83  77
4 04.01.1961  980.4  982.6  986.6  983.2  5.5  1.5 -1.3 3.4 5.4 2.6  84  85  92
5 05.01.1961  993.6 1001.4 1003.9  999.6  6.4  1.7 -0.4 4.4 5.9 2.0  77  65  82
  rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd  sd          vr1          vr2
1  91  C  C  C   0   0   0 0.0   8   8   8 8.0 0.0 0.5 bis 1 km   2 bis 4 km
2  89  S  S NO   2   1   2 1.7   7   8   8 7.7 0.0 10 bis 20 km  4 bis 10 km
3  79  S  S SO   5   1   1 2.3   8   7   7 7.3 0.0  4 bis 10 km 20 bis 50 km
4  87  S  S  S   1   2   2 1.7   7   6   7 6.7 0.3  4 bis 10 km 20 bis 50 km
5  75 SW SW  S   3   3   2 2.7   6   4   7 5.7 2.5  4 bis 10 km 20 bis 50 km
           vr3 pr1  prt1 pr2  prt2 pr3  prt3 pd.1   pdt sh wsdmax
1   1 bis 2 km 0.0     - 0.0     - 0.0     -  0.0     -  0    3.5
2   2 bis 4 km 0.0     - 0.7 Regen 1.7 Regen  7.9 Regen  0    6.9
3 10 bis 20 km 5.5 Regen 1.4 Regen 0.0     -  1.4 Regen  0   18.4
4  4 bis 10 km 0.0 Regen 0.6 Regen 4.3 Regen  5.0 Regen  0   14.1
5  4 bis 10 km 0.1 Regen 0.0 Regen 0.0     -  0.0 Regen  0   15.0

7.2.3 Code 7486fb94

load("data/c55d0a7b.RData")
year=format(d$da1, format="%Y")
year[1:5]
d=subset(d,year=="1990")
d[1:5,]
rownames(d)=1:nrow(d)
save(d,file="data/7486fb94.RData")

load("data/c55d0a7b.RData")
year=format(d$da1, format="%Y")
year[1:5]

[1] "1961" "1961" "1961" "1961" "1961"

d=subset(d,year=="1990")
d[1:5,]

              da     p1     p2     p3     pd tmax tmin emin   t1   t2   t3 rm1
10593 01.01.1990 1005.4 1004.7 1005.7 1005.3 -1.0 -2.0 -1.6 -1.3 -1.7 -2.0  98
10594 02.01.1990 1006.9 1007.8 1009.0 1007.9 -1.0 -2.4 -2.5 -2.4 -1.9 -1.0  96
10595 03.01.1990 1008.5 1008.2 1009.7 1008.8  0.3 -1.0 -1.6  0.2 -0.2  0.3  95
10596 04.01.1990 1010.1 1010.1 1010.9 1010.4  0.7 -0.2 -0.6  0.1  0.7  0.0  73
10597 05.01.1990 1010.2 1010.7 1012.3 1011.1  1.5 -0.4 -0.8  0.2  1.2  1.3  91
      rm2 rm3 rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd sd         vr1
10593  94  94  95  O  O SO   1   1   1 1.0   8   8   8   8  0  2 bis 4 km
10594  94  96  95  S  S SO   1   1   1 1.0   8   8   8   8  0  2 bis 4 km
10595  84  72  84  O  O  O   2   2   2 2.0   8   8   8   8  0 4 bis 10 km
10596  74  82  76 SO  S  S   2   2   2 2.0   8   8   8   8  0 4 bis 10 km
10597  95  95  94  S  C NO   1   0   1 0.7   8   8   8   8  0 4 bis 10 km
              vr2         vr3 pr1 prt1 pr2   prt2 pr3 prt3 pd.1    pdt sh
10593  2 bis 4 km  2 bis 4 km   0    -   0      -   0    -    0      -  0
10594  1 bis 2 km  2 bis 4 km   0    -   0      -   0    -    0      -  0
10595 4 bis 10 km 4 bis 10 km   0    -   0 Schnee   0    -    0 Schnee  0
10596 4 bis 10 km 4 bis 10 km   0    -   0      -   0    -    0      -  0
10597  2 bis 4 km 4 bis 10 km   0    -   0      -   0    -    0      -  0
      wsdmax        da1
10593    2.9 1990-01-01
10594    2.0 1990-01-02
10595    5.8 1990-01-03
10596    4.4 1990-01-04
10597    3.0 1990-01-05

rownames(d)=1:nrow(d)
save(d,file="data/7486fb94.RData")

7.3 Knowledge database

7.3.1 Data

7.3.1.1 Code dc602b67

load("data/7486fb94.RData")
DT::datatable(
  d   # USER INPUT: replace d below by the appropriate data frame 
  , extensions = 'FixedColumns'
  ,options = list(
    dom = 't',
    scrollX = TRUE,
    scrollY = 400,
    fixedColumns = TRUE
    ,pageLength = nrow(d)
    # ,lengthMenu = c( 7,20)
  )
)

7.3.1.2 Code 16879b5a

load("data/7486fb94.RData")
str(d)
class(d$da)
class(d$da1)
class(d$t1)
class(d$ws1)
class(d$w1)
d$w1f=as.factor(d$w1)
class(d$w1f)
levels(d$w1f)
class(d$rm1)
d$rm1n=as.numeric(d$rm1)
class(d$rm1n)
t2c=ggplot2::cut_interval(d$t2,n = 3)
class(t2c)
table(t2c)
data.frame(t2=d$t2,t2c=t2c)
t2c=ggplot2::cut_number(d$t2,n = 3)
table(t2c)
summary(d)

load("data/7486fb94.RData")
str(d)

'data.frame':   365 obs. of  41 variables:
 $ da    : chr  "01.01.1990" "02.01.1990" "03.01.1990" "04.01.1990" ...
 $ p1    : num  1005 1007 1008 1010 1010 ...
 $ p2    : num  1005 1008 1008 1010 1011 ...
 $ p3    : num  1006 1009 1010 1011 1012 ...
 $ pd    : num  1005 1008 1009 1010 1011 ...
 $ tmax  : num  -1 -1 0.3 0.7 1.5 1.3 1.7 1.1 3.3 3 ...
 $ tmin  : num  -2 -2.4 -1 -0.2 -0.4 0.2 1 -1 1 1.6 ...
 $ emin  : num  -1.6 -2.5 -1.6 -0.6 -0.8 0.5 0.3 -1.3 0.6 1.3 ...
 $ t1    : num  -1.3 -2.4 0.2 0.1 0.2 0.8 1.6 -0.8 1.8 1.8 ...
 $ t2    : num  -1.7 -1.9 -0.2 0.7 1.2 0.5 1.7 0.9 3.3 2.1 ...
 $ t3    : num  -2 -1 0.3 0 1.3 1 1 1.1 3 2.6 ...
 $ rm1   : int  98 96 95 73 91 91 80 96 97 93 ...
 $ rm2   : int  94 94 84 74 95 88 83 91 94 88 ...
 $ rm3   : int  94 96 72 82 95 83 86 93 94 89 ...
 $ rmd   : int  95 95 84 76 94 87 83 93 95 90 ...
 $ w1    : chr  "O" "S" "O" "SO" ...
 $ w2    : chr  "O" "S" "O" "S" ...
 $ w3    : chr  "SO" "SO" "O" "S" ...
 $ ws1   : int  1 1 2 2 1 2 2 2 1 2 ...
 $ ws2   : int  1 1 2 2 0 2 2 1 0 1 ...
 $ ws3   : int  1 1 2 2 1 2 2 1 1 1 ...
 $ wsd   : num  1 1 2 2 0.7 2 2 1.3 0.7 1.3 ...
 $ cc1   : int  8 8 8 8 8 8 8 8 8 8 ...
 $ cc2   : int  8 8 8 8 8 8 8 8 8 8 ...
 $ cc3   : int  8 8 8 8 8 8 8 8 8 8 ...
 $ ccd   : num  8 8 8 8 8 8 8 8 8 8 ...
 $ sd    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ vr1   : chr  "2 bis 4 km" "2 bis 4 km" "4 bis 10 km" "4 bis 10 km" ...
 $ vr2   : chr  "2 bis 4 km" "1 bis 2 km" "4 bis 10 km" "4 bis 10 km" ...
 $ vr3   : chr  "2 bis 4 km" "2 bis 4 km" "4 bis 10 km" "4 bis 10 km" ...
 $ pr1   : num  0 0 0 0 0 0 0 0 0.5 0 ...
 $ prt1  : chr  "-" "-" "-" "-" ...
 $ pr2   : num  0 0 0 0 0 0 0 0 0.3 0 ...
 $ prt2  : chr  "-" "-" "Schnee" "-" ...
 $ pr3   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ prt3  : chr  "-" "-" "-" "-" ...
 $ pd.1  : num  0 0 0 0 0 0 0 0.5 0.3 0 ...
 $ pdt   : chr  "-" "-" "Schnee" "-" ...
 $ sh    : int  0 0 0 0 0 0 0 0 0 0 ...
 $ wsdmax: num  2.9 2 5.8 4.4 3 4.9 3.4 2.1 2.1 3.7 ...
 $ da1   : Date, format: "1990-01-01" "1990-01-02" ...

class(d$da)

[1] "character"

class(d$da1)

[1] "Date"

class(d$t1)

[1] "numeric"

class(d$ws1)

[1] "integer"

class(d$w1)

[1] "character"

d$w1f=as.factor(d$w1)
class(d$w1f)

[1] "factor"

levels(d$w1f)

[1] "C"  "N"  "NO" "NW" "O"  "S"  "SO" "SW" "W"

class(d$rm1)

[1] "integer"

d$rm1n=as.numeric(d$rm1)
class(d$rm1n)

[1] "numeric"

t2c=ggplot2::cut_interval(d$t2,n = 3)
class(t2c)

[1] "factor"

table(t2c)

t2c
 [-2.1,9.7]  (9.7,21.5] (21.5,33.3] 
        108         189          68

data.frame(t2=d$t2,t2c=t2c)[16:18,]

     t2        t2c
16  7.8 [-2.1,9.7]
17 10.3 (9.7,21.5]
18  6.6 [-2.1,9.7]

t2c=ggplot2::cut_number(d$t2,n = 3)
table(t2c)

t2c
[-2.1,10.4] (10.4,17.8] (17.8,33.3] 
        122         122         121

summary(d)

      da                  p1               p2               p3        
 Length:365         Min.   : 971.0   Min.   : 974.3   Min.   : 969.0  
 Class :character   1st Qu.: 998.6   1st Qu.: 997.8   1st Qu.: 997.5  
 Mode  :character   Median :1003.7   Median :1003.3   Median :1003.8  
                    Mean   :1003.4   Mean   :1002.7   Mean   :1003.2  
                    3rd Qu.:1010.1   3rd Qu.:1009.0   3rd Qu.:1009.4  
                    Max.   :1026.6   Max.   :1026.9   Max.   :1028.3  
                                                                      
       pd              tmax            tmin             emin       
 Min.   : 971.4   Min.   :-1.00   Min.   :-5.900   Min.   :-8.500  
 1st Qu.: 997.9   1st Qu.: 9.40   1st Qu.: 2.500   1st Qu.: 0.500  
 Median :1003.8   Median :15.70   Median : 6.800   Median : 5.000  
 Mean   :1003.1   Mean   :15.34   Mean   : 6.941   Mean   : 4.998  
 3rd Qu.:1009.2   3rd Qu.:20.80   3rd Qu.:11.000   3rd Qu.: 9.200  
 Max.   :1026.6   Max.   :33.90   Max.   :20.400   Max.   :18.200  
                                                                   
       t1               t2              t3             rm1        
 Min.   :-5.600   Min.   :-2.10   Min.   :-3.20   Min.   : 40.00  
 1st Qu.: 3.800   1st Qu.: 8.10   1st Qu.: 5.50   1st Qu.: 75.00  
 Median : 8.800   Median :14.80   Median :10.40   Median : 84.00  
 Mean   : 8.723   Mean   :14.31   Mean   :10.48   Mean   : 83.13  
 3rd Qu.:13.400   3rd Qu.:19.70   3rd Qu.:15.00   3rd Qu.: 93.00  
 Max.   :22.900   Max.   :33.30   Max.   :26.10   Max.   :100.00  
                                                                  
      rm2             rm3            rmd            w1           
 Min.   :18.00   Min.   :29.0   Min.   :33.0   Length:365        
 1st Qu.:44.00   1st Qu.:65.0   1st Qu.:63.0   Class :character  
 Median :57.00   Median :80.0   Median :74.0   Mode  :character  
 Mean   :58.68   Mean   :75.6   Mean   :72.5                     
 3rd Qu.:72.00   3rd Qu.:90.0   3rd Qu.:82.0                     
 Max.   :97.00   Max.   :98.0   Max.   :97.0                     
                                                                 
      w2                 w3                 ws1             ws2       
 Length:365         Length:365         Min.   :0.000   Min.   :0.000  
 Class :character   Class :character   1st Qu.:1.000   1st Qu.:1.000  
 Mode  :character   Mode  :character   Median :1.000   Median :2.000  
                                       Mean   :1.499   Mean   :2.381  
                                       3rd Qu.:2.000   3rd Qu.:3.000  
                                       Max.   :5.000   Max.   :6.000  
                                                                      
      ws3             wsd             cc1             cc2       
 Min.   :0.000   Min.   :0.300   Min.   :0.000   Min.   :0.000  
 1st Qu.:1.000   1st Qu.:1.300   1st Qu.:3.000   1st Qu.:4.000  
 Median :2.000   Median :1.700   Median :6.000   Median :7.000  
 Mean   :1.808   Mean   :1.894   Mean   :5.263   Mean   :5.493  
 3rd Qu.:2.000   3rd Qu.:2.300   3rd Qu.:8.000   3rd Qu.:8.000  
 Max.   :5.000   Max.   :5.000   Max.   :8.000   Max.   :8.000  
                                                                
      cc3             ccd              sd             vr1           
 Min.   :0.000   Min.   :0.000   Min.   : 0.000   Length:365        
 1st Qu.:1.000   1st Qu.:3.300   1st Qu.: 0.400   Class :character  
 Median :6.000   Median :5.700   Median : 4.300   Mode  :character  
 Mean   :4.564   Mean   :5.107   Mean   : 4.915                     
 3rd Qu.:8.000   3rd Qu.:7.000   3rd Qu.: 8.100                     
 Max.   :8.000   Max.   :8.000   Max.   :15.200                     
                                                                    
     vr2                vr3                 pr1              prt1          
 Length:365         Length:365         Min.   : 0.0000   Length:365        
 Class :character   Class :character   1st Qu.: 0.0000   Class :character  
 Mode  :character   Mode  :character   Median : 0.0000   Mode  :character  
                                       Mean   : 0.4277                     
                                       3rd Qu.: 0.1000                     
                                       Max.   :15.8000                     
                                                                           
      pr2              prt2                pr3              prt3          
 Min.   : 0.0000   Length:365         Min.   : 0.0000   Length:365        
 1st Qu.: 0.0000   Class :character   1st Qu.: 0.0000   Class :character  
 Median : 0.0000   Mode  :character   Median : 0.0000   Mode  :character  
 Mean   : 0.4625                      Mean   : 0.4759                     
 3rd Qu.: 0.0000                      3rd Qu.: 0.0000                     
 Max.   :18.1000                      Max.   :19.6000                     
                                                                          
      pd.1            pdt                  sh             wsdmax      
 Min.   : 0.000   Length:365         Min.   :0.0000   Min.   : 1.300  
 1st Qu.: 0.000   Class :character   1st Qu.:0.0000   1st Qu.: 5.100  
 Median : 0.000   Mode  :character   Median :0.0000   Median : 7.300  
 Mean   : 1.366                      Mean   :0.1041   Mean   : 8.071  
 3rd Qu.: 1.100                      3rd Qu.:0.0000   3rd Qu.:10.000  
 Max.   :24.900                      Max.   :9.0000   Max.   :28.400  
                                                                      
      da1                  w1f           rm1n       
 Min.   :1990-01-01   S      :108   Min.   : 40.00  
 1st Qu.:1990-04-02   W      : 53   1st Qu.: 75.00  
 Median :1990-07-02   SW     : 46   Median : 84.00  
 Mean   :1990-07-02   NW     : 39   Mean   : 83.13  
 3rd Qu.:1990-10-01   C      : 38   3rd Qu.: 93.00  
 Max.   :1990-12-31   NO     : 30   Max.   :100.00  
                      (Other): 51

7.3.1.3 Code a8162e80

load("data/7486fb94.RData")
d=d[1:5,]
d
d$t1
d[,]
d[,c("t1")]
d[,c("t1","t2","t3")]
d[,c(9,10,11)]
d[,9:11]
d[2,9:11]
d[2:3,9:11]
d[c(1,3,5),9:11]
d[1:2,]
d[c(F,T,F,T,F),]
cols=c(rep(F,8),c(T,T,T),rep(F,ncol(d)-11))
cols
d[,cols]
d[1]
d[2]

load("data/7486fb94.RData")
d=d[1:5,]
d

          da     p1     p2     p3     pd tmax tmin emin   t1   t2   t3 rm1 rm2
1 01.01.1990 1005.4 1004.7 1005.7 1005.3 -1.0 -2.0 -1.6 -1.3 -1.7 -2.0  98  94
2 02.01.1990 1006.9 1007.8 1009.0 1007.9 -1.0 -2.4 -2.5 -2.4 -1.9 -1.0  96  94
3 03.01.1990 1008.5 1008.2 1009.7 1008.8  0.3 -1.0 -1.6  0.2 -0.2  0.3  95  84
4 04.01.1990 1010.1 1010.1 1010.9 1010.4  0.7 -0.2 -0.6  0.1  0.7  0.0  73  74
5 05.01.1990 1010.2 1010.7 1012.3 1011.1  1.5 -0.4 -0.8  0.2  1.2  1.3  91  95
  rm3 rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd sd         vr1         vr2
1  94  95  O  O SO   1   1   1 1.0   8   8   8   8  0  2 bis 4 km  2 bis 4 km
2  96  95  S  S SO   1   1   1 1.0   8   8   8   8  0  2 bis 4 km  1 bis 2 km
3  72  84  O  O  O   2   2   2 2.0   8   8   8   8  0 4 bis 10 km 4 bis 10 km
4  82  76 SO  S  S   2   2   2 2.0   8   8   8   8  0 4 bis 10 km 4 bis 10 km
5  95  94  S  C NO   1   0   1 0.7   8   8   8   8  0 4 bis 10 km  2 bis 4 km
          vr3 pr1 prt1 pr2   prt2 pr3 prt3 pd.1    pdt sh wsdmax        da1
1  2 bis 4 km   0    -   0      -   0    -    0      -  0    2.9 1990-01-01
2  2 bis 4 km   0    -   0      -   0    -    0      -  0    2.0 1990-01-02
3 4 bis 10 km   0    -   0 Schnee   0    -    0 Schnee  0    5.8 1990-01-03
4 4 bis 10 km   0    -   0      -   0    -    0      -  0    4.4 1990-01-04
5 4 bis 10 km   0    -   0      -   0    -    0      -  0    3.0 1990-01-05

d$t1

[1] -1.3 -2.4  0.2  0.1  0.2

d[,]

          da     p1     p2     p3     pd tmax tmin emin   t1   t2   t3 rm1 rm2
1 01.01.1990 1005.4 1004.7 1005.7 1005.3 -1.0 -2.0 -1.6 -1.3 -1.7 -2.0  98  94
2 02.01.1990 1006.9 1007.8 1009.0 1007.9 -1.0 -2.4 -2.5 -2.4 -1.9 -1.0  96  94
3 03.01.1990 1008.5 1008.2 1009.7 1008.8  0.3 -1.0 -1.6  0.2 -0.2  0.3  95  84
4 04.01.1990 1010.1 1010.1 1010.9 1010.4  0.7 -0.2 -0.6  0.1  0.7  0.0  73  74
5 05.01.1990 1010.2 1010.7 1012.3 1011.1  1.5 -0.4 -0.8  0.2  1.2  1.3  91  95
  rm3 rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd sd         vr1         vr2
1  94  95  O  O SO   1   1   1 1.0   8   8   8   8  0  2 bis 4 km  2 bis 4 km
2  96  95  S  S SO   1   1   1 1.0   8   8   8   8  0  2 bis 4 km  1 bis 2 km
3  72  84  O  O  O   2   2   2 2.0   8   8   8   8  0 4 bis 10 km 4 bis 10 km
4  82  76 SO  S  S   2   2   2 2.0   8   8   8   8  0 4 bis 10 km 4 bis 10 km
5  95  94  S  C NO   1   0   1 0.7   8   8   8   8  0 4 bis 10 km  2 bis 4 km
          vr3 pr1 prt1 pr2   prt2 pr3 prt3 pd.1    pdt sh wsdmax        da1
1  2 bis 4 km   0    -   0      -   0    -    0      -  0    2.9 1990-01-01
2  2 bis 4 km   0    -   0      -   0    -    0      -  0    2.0 1990-01-02
3 4 bis 10 km   0    -   0 Schnee   0    -    0 Schnee  0    5.8 1990-01-03
4 4 bis 10 km   0    -   0      -   0    -    0      -  0    4.4 1990-01-04
5 4 bis 10 km   0    -   0      -   0    -    0      -  0    3.0 1990-01-05

d[,c("t1")]

[1] -1.3 -2.4  0.2  0.1  0.2

d[,c("t1","t2","t3")]

    t1   t2   t3
1 -1.3 -1.7 -2.0
2 -2.4 -1.9 -1.0
3  0.2 -0.2  0.3
4  0.1  0.7  0.0
5  0.2  1.2  1.3

d[,c(9,10,11)]

    t1   t2   t3
1 -1.3 -1.7 -2.0
2 -2.4 -1.9 -1.0
3  0.2 -0.2  0.3
4  0.1  0.7  0.0
5  0.2  1.2  1.3

d[,9:11]

    t1   t2   t3
1 -1.3 -1.7 -2.0
2 -2.4 -1.9 -1.0
3  0.2 -0.2  0.3
4  0.1  0.7  0.0
5  0.2  1.2  1.3

d[2,9:11]

    t1   t2 t3
2 -2.4 -1.9 -1

d[2:3,9:11]

    t1   t2   t3
2 -2.4 -1.9 -1.0
3  0.2 -0.2  0.3

d[c(1,3,5),9:11]

    t1   t2   t3
1 -1.3 -1.7 -2.0
3  0.2 -0.2  0.3
5  0.2  1.2  1.3

d[1:2,]

          da     p1     p2     p3     pd tmax tmin emin   t1   t2 t3 rm1 rm2
1 01.01.1990 1005.4 1004.7 1005.7 1005.3   -1 -2.0 -1.6 -1.3 -1.7 -2  98  94
2 02.01.1990 1006.9 1007.8 1009.0 1007.9   -1 -2.4 -2.5 -2.4 -1.9 -1  96  94
  rm3 rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd sd        vr1        vr2
1  94  95  O  O SO   1   1   1   1   8   8   8   8  0 2 bis 4 km 2 bis 4 km
2  96  95  S  S SO   1   1   1   1   8   8   8   8  0 2 bis 4 km 1 bis 2 km
         vr3 pr1 prt1 pr2 prt2 pr3 prt3 pd.1 pdt sh wsdmax        da1
1 2 bis 4 km   0    -   0    -   0    -    0   -  0    2.9 1990-01-01
2 2 bis 4 km   0    -   0    -   0    -    0   -  0    2.0 1990-01-02

d[c(F,T,F,T,F),]

          da     p1     p2     p3     pd tmax tmin emin   t1   t2 t3 rm1 rm2
2 02.01.1990 1006.9 1007.8 1009.0 1007.9 -1.0 -2.4 -2.5 -2.4 -1.9 -1  96  94
4 04.01.1990 1010.1 1010.1 1010.9 1010.4  0.7 -0.2 -0.6  0.1  0.7  0  73  74
  rm3 rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd sd         vr1         vr2
2  96  95  S  S SO   1   1   1   1   8   8   8   8  0  2 bis 4 km  1 bis 2 km
4  82  76 SO  S  S   2   2   2   2   8   8   8   8  0 4 bis 10 km 4 bis 10 km
          vr3 pr1 prt1 pr2 prt2 pr3 prt3 pd.1 pdt sh wsdmax        da1
2  2 bis 4 km   0    -   0    -   0    -    0   -  0    2.0 1990-01-02
4 4 bis 10 km   0    -   0    -   0    -    0   -  0    4.4 1990-01-04

cols=c(rep(F,8),c(T,T,T),rep(F,ncol(d)-11))
cols

 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE FALSE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[37] FALSE FALSE FALSE FALSE FALSE

d[,cols]

    t1   t2   t3
1 -1.3 -1.7 -2.0
2 -2.4 -1.9 -1.0
3  0.2 -0.2  0.3
4  0.1  0.7  0.0
5  0.2  1.2  1.3

d[1]

          da
1 01.01.1990
2 02.01.1990
3 03.01.1990
4 04.01.1990
5 05.01.1990

d[2]

7.3.1.4 Code 854728e5

load("data/7486fb94.RData")
d=d[151:155,c("t1","t2","w1","w2","w3")]
d
d$w2==d$w3
data.frame(d,l1=d$w2==d$w3)
d[d$w2==d$w3,]
subset(d,d$w2==d$w3)
!d$w2==d$w3
data.frame(d,l2=!d$w2==d$w3)
!d$w2==d$w3 & d$t1>13
data.frame(d,l3=!d$w2==d$w3 & d$t1>13)
!d$w2==d$w3 | d$t1>13
data.frame(d,l4=!d$w2==d$w3 | d$t1>13)
colnames(d)
colnames(d) %in% c("w1","w3")
d[,colnames(d) %in% c("w1","w3")]
d[,!(colnames(d) %in% c("w1","w3"))]

load("data/7486fb94.RData")
d=d[151:155,c("t1","t2","w1","w2","w3")]
d

      t1   t2 w1 w2 w3
151 10.6 24.0  C SO NW
152 13.4 26.2 NW  S  W
153 14.8 19.9  S SW SW
154 12.8 16.7 SO  S  S
155 13.1 16.0 SW  W SW

d$w2==d$w3

[1] FALSE FALSE  TRUE  TRUE FALSE

data.frame(d,l1=d$w2==d$w3)

      t1   t2 w1 w2 w3    l1
151 10.6 24.0  C SO NW FALSE
152 13.4 26.2 NW  S  W FALSE
153 14.8 19.9  S SW SW  TRUE
154 12.8 16.7 SO  S  S  TRUE
155 13.1 16.0 SW  W SW FALSE

d[d$w2==d$w3,]

      t1   t2 w1 w2 w3
153 14.8 19.9  S SW SW
154 12.8 16.7 SO  S  S

subset(d,d$w2==d$w3)

      t1   t2 w1 w2 w3
153 14.8 19.9  S SW SW
154 12.8 16.7 SO  S  S

!d$w2==d$w3

[1]  TRUE  TRUE FALSE FALSE  TRUE

data.frame(d,l2=!d$w2==d$w3)

      t1   t2 w1 w2 w3    l2
151 10.6 24.0  C SO NW  TRUE
152 13.4 26.2 NW  S  W  TRUE
153 14.8 19.9  S SW SW FALSE
154 12.8 16.7 SO  S  S FALSE
155 13.1 16.0 SW  W SW  TRUE

!d$w2==d$w3 & d$t1>13

[1] FALSE  TRUE FALSE FALSE  TRUE

data.frame(d,l3=!d$w2==d$w3 & d$t1>13)

      t1   t2 w1 w2 w3    l3
151 10.6 24.0  C SO NW FALSE
152 13.4 26.2 NW  S  W  TRUE
153 14.8 19.9  S SW SW FALSE
154 12.8 16.7 SO  S  S FALSE
155 13.1 16.0 SW  W SW  TRUE

!d$w2==d$w3 | d$t1>13

[1]  TRUE  TRUE  TRUE FALSE  TRUE

data.frame(d,l4=!d$w2==d$w3 | d$t1>13)

      t1   t2 w1 w2 w3    l4
151 10.6 24.0  C SO NW  TRUE
152 13.4 26.2 NW  S  W  TRUE
153 14.8 19.9  S SW SW  TRUE
154 12.8 16.7 SO  S  S FALSE
155 13.1 16.0 SW  W SW  TRUE

colnames(d)

[1] "t1" "t2" "w1" "w2" "w3"

colnames(d) %in% c("w1","w3")

[1] FALSE FALSE  TRUE FALSE  TRUE

d[,colnames(d) %in% c("w1","w3")]

    w1 w3
151  C NW
152 NW  W
153  S SW
154 SO  S
155 SW SW

d[,!(colnames(d) %in% c("w1","w3"))]

      t1   t2 w2
151 10.6 24.0 SO
152 13.4 26.2  S
153 14.8 19.9 SW
154 12.8 16.7  S
155 13.1 16.0  W

7.3.1.5 Code 34d67994

d=read.csv2("data/weather2.csv")
d[1:5,]
class(d$da)
d$da=as.Date(d$da,"%d.%m.%Y")
class(d$da)
d[1:5,]
data.frame(Date=d$da,Day1=format(d$da, format="%d"))[c(1:3,30:34),]
data.frame(Date=d$da,Day2=format(d$da, format="%j"))[c(1:3,364:369),]
data.frame(Date=d$da,Day3=weekdays(d$da))[1:3,]
data.frame(Date=d$da,Day4=format(d$da, format="%u"))[1:9,]
# %U: Week of the year as decimal number (00–53) using Sunday as the first day 1 of the week (and typically with the first Sunday of the year as day 1 of week 1). The US convention.
data.frame(Date=d$da,Week1=format(d$da, format="%U"))[c(1:3,32:34,364:372),]
# %W: Week of the year as decimal number (00–53) using Monday as the first day of week (and typically with the first Monday of the year as day 1 of week 1). The UK convention.
data.frame(Date=d$da,Week2=format(d$da, format="%W"))[c(1:3,32:34,364:369),]
data.frame(Date=d$da,Month1=months(d$da))[c(1:3,32:34),]
data.frame(Date=d$da,Month2=format(d$da, format="%m"))[c(1:3,32:34),]
data.frame(Date=d$da,Quarter=quarters(d$da))[c(1:3,91:93),]
data.frame(Date=d$da,Years=format(d$da, format="%Y"))[c(1:3,366:368),]

Help pages: as.Date, strptime

d=read.csv2("data/weather2.csv")
d[1:5,]

          da     p1     p2     p3     pd tmax tmin emin  t1  t2  t3 rm1 rm2 rm3
1 01.01.1961 1001.4 1001.5 1001.8 1001.6  2.8  1.4  0.0 1.5 2.6 2.4  95  90  87
2 02.01.1961  996.7  992.0  985.8  991.5  4.1  0.9 -1.1 3.1 2.6 3.7  77  94  95
3 03.01.1961  973.5  980.6  981.4  978.5  8.6  2.7  3.0 8.4 4.9 3.5  78  83  77
4 04.01.1961  980.4  982.6  986.6  983.2  5.5  1.5 -1.3 3.4 5.4 2.6  84  85  92
5 05.01.1961  993.6 1001.4 1003.9  999.6  6.4  1.7 -0.4 4.4 5.9 2.0  77  65  82
  rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd  sd          vr1          vr2
1  91  C  C  C   0   0   0 0.0   8   8   8 8.0 0.0 0.5 bis 1 km   2 bis 4 km
2  89  S  S NO   2   1   2 1.7   7   8   8 7.7 0.0 10 bis 20 km  4 bis 10 km
3  79  S  S SO   5   1   1 2.3   8   7   7 7.3 0.0  4 bis 10 km 20 bis 50 km
4  87  S  S  S   1   2   2 1.7   7   6   7 6.7 0.3  4 bis 10 km 20 bis 50 km
5  75 SW SW  S   3   3   2 2.7   6   4   7 5.7 2.5  4 bis 10 km 20 bis 50 km
           vr3 pr1  prt1 pr2  prt2 pr3  prt3 pd.1   pdt sh wsdmax
1   1 bis 2 km 0.0     - 0.0     - 0.0     -  0.0     -  0    3.5
2   2 bis 4 km 0.0     - 0.7 Regen 1.7 Regen  7.9 Regen  0    6.9
3 10 bis 20 km 5.5 Regen 1.4 Regen 0.0     -  1.4 Regen  0   18.4
4  4 bis 10 km 0.0 Regen 0.6 Regen 4.3 Regen  5.0 Regen  0   14.1
5  4 bis 10 km 0.1 Regen 0.0 Regen 0.0     -  0.0 Regen  0   15.0

class(d$da)

[1] "character"

d$da=as.Date(d$da,"%d.%m.%Y")
class(d$da)

[1] "Date"

d[1:5,]

          da     p1     p2     p3     pd tmax tmin emin  t1  t2  t3 rm1 rm2 rm3
1 1961-01-01 1001.4 1001.5 1001.8 1001.6  2.8  1.4  0.0 1.5 2.6 2.4  95  90  87
2 1961-01-02  996.7  992.0  985.8  991.5  4.1  0.9 -1.1 3.1 2.6 3.7  77  94  95
3 1961-01-03  973.5  980.6  981.4  978.5  8.6  2.7  3.0 8.4 4.9 3.5  78  83  77
4 1961-01-04  980.4  982.6  986.6  983.2  5.5  1.5 -1.3 3.4 5.4 2.6  84  85  92
5 1961-01-05  993.6 1001.4 1003.9  999.6  6.4  1.7 -0.4 4.4 5.9 2.0  77  65  82
  rmd w1 w2 w3 ws1 ws2 ws3 wsd cc1 cc2 cc3 ccd  sd          vr1          vr2
1  91  C  C  C   0   0   0 0.0   8   8   8 8.0 0.0 0.5 bis 1 km   2 bis 4 km
2  89  S  S NO   2   1   2 1.7   7   8   8 7.7 0.0 10 bis 20 km  4 bis 10 km
3  79  S  S SO   5   1   1 2.3   8   7   7 7.3 0.0  4 bis 10 km 20 bis 50 km
4  87  S  S  S   1   2   2 1.7   7   6   7 6.7 0.3  4 bis 10 km 20 bis 50 km
5  75 SW SW  S   3   3   2 2.7   6   4   7 5.7 2.5  4 bis 10 km 20 bis 50 km
           vr3 pr1  prt1 pr2  prt2 pr3  prt3 pd.1   pdt sh wsdmax
1   1 bis 2 km 0.0     - 0.0     - 0.0     -  0.0     -  0    3.5
2   2 bis 4 km 0.0     - 0.7 Regen 1.7 Regen  7.9 Regen  0    6.9
3 10 bis 20 km 5.5 Regen 1.4 Regen 0.0     -  1.4 Regen  0   18.4
4  4 bis 10 km 0.0 Regen 0.6 Regen 4.3 Regen  5.0 Regen  0   14.1
5  4 bis 10 km 0.1 Regen 0.0 Regen 0.0     -  0.0 Regen  0   15.0

data.frame(Date=d$da,Day1=format(d$da, format="%d"))[c(1:3,30:34),]

         Date Day1
1  1961-01-01   01
2  1961-01-02   02
3  1961-01-03   03
30 1961-01-30   30
31 1961-01-31   31
32 1961-02-01   01
33 1961-02-02   02
34 1961-02-03   03

data.frame(Date=d$da,Day2=format(d$da, format="%j"))[c(1:3,364:369),]

          Date Day2
1   1961-01-01  001
2   1961-01-02  002
3   1961-01-03  003
364 1961-12-30  364
365 1961-12-31  365
366 1962-01-01  001
367 1962-01-02  002
368 1962-01-03  003
369 1962-01-04  004

data.frame(Date=d$da,Day3=weekdays(d$da))[1:3,]

        Date     Day3
1 1961-01-01  Sonntag
2 1961-01-02   Montag
3 1961-01-03 Dienstag

data.frame(Date=d$da,Day4=format(d$da, format="%u"))[1:9,]

        Date Day4
1 1961-01-01    7
2 1961-01-02    1
3 1961-01-03    2
4 1961-01-04    3
5 1961-01-05    4
6 1961-01-06    5
7 1961-01-07    6
8 1961-01-08    7
9 1961-01-09    1

# %U: Week of the year as decimal number (00–53) using Sunday as the first day 1 of the week (and typically with the first Sunday of the year as day 1 of week 1). The US convention.
data.frame(Date=d$da,Week1=format(d$da, format="%U"))[c(1:3,32:34,364:372),]

          Date Week1
1   1961-01-01    01
2   1961-01-02    01
3   1961-01-03    01
32  1961-02-01    05
33  1961-02-02    05
34  1961-02-03    05
364 1961-12-30    52
365 1961-12-31    53
366 1962-01-01    00
367 1962-01-02    00
368 1962-01-03    00
369 1962-01-04    00
370 1962-01-05    00
371 1962-01-06    00
372 1962-01-07    01

# %W: Week of the year as decimal number (00–53) using Monday as the first day of week (and typically with the first Monday of the year as day 1 of week 1). The UK convention.
data.frame(Date=d$da,Week2=format(d$da, format="%W"))[c(1:3,32:34,364:369),]

          Date Week2
1   1961-01-01    00
2   1961-01-02    01
3   1961-01-03    01
32  1961-02-01    05
33  1961-02-02    05
34  1961-02-03    05
364 1961-12-30    52
365 1961-12-31    52
366 1962-01-01    01
367 1962-01-02    01
368 1962-01-03    01
369 1962-01-04    01

data.frame(Date=d$da,Month1=months(d$da))[c(1:3,32:34),]

         Date  Month1
1  1961-01-01  Januar
2  1961-01-02  Januar
3  1961-01-03  Januar
32 1961-02-01 Februar
33 1961-02-02 Februar
34 1961-02-03 Februar

data.frame(Date=d$da,Month2=format(d$da, format="%m"))[c(1:3,32:34),]

         Date Month2
1  1961-01-01     01
2  1961-01-02     01
3  1961-01-03     01
32 1961-02-01     02
33 1961-02-02     02
34 1961-02-03     02

data.frame(Date=d$da,Quarter=quarters(d$da))[c(1:3,91:93),]

         Date Quarter
1  1961-01-01      Q1
2  1961-01-02      Q1
3  1961-01-03      Q1
91 1961-04-01      Q2
92 1961-04-02      Q2
93 1961-04-03      Q2

data.frame(Date=d$da,Years=format(d$da, format="%Y"))[c(1:3,366:368),]

          Date Years
1   1961-01-01  1961
2   1961-01-02  1961
3   1961-01-03  1961
366 1962-01-01  1962
367 1962-01-02  1962
368 1962-01-03  1962

7.3.2 Plots

7.3.2.1 Code 002a7650

load("data/7486fb94.RData")
table(d$w1)
barplot(table(d$w1))

load("data/7486fb94.RData")
table(d$w1)


  C   N  NO  NW   O   S  SO  SW   W 
 38  14  30  39  25 108  12  46  53

barplot(table(d$w1))

7.3.2.2 Code ae53a438

load("data/7486fb94.RData")
library(ggplot2)
p=ggplot(d,aes(w1))
p=p+geom_bar()
p=p+xlab("Wind direction 7.30 am")
p=p+ylab("Frequency")
p

7.3.2.3 Code 855a5ed8

load("data/7486fb94.RData")
library(ggplot2)
p=ggplot(d,aes(w1,t1))
p=p+xlab("Wind direction 7.30 am")
p=p+ylab("Temperature 7.30 am")
p=p+geom_boxplot(fill="red")
p
plotly::ggplotly(p)

load("data/7486fb94.RData")
library(ggplot2)
p=ggplot(d,aes(w1,t1))
p=p+xlab("Wind direction 7.30 am")
p=p+ylab("Temperature 7.30 am")
p=p+geom_boxplot(fill="red")
p

plotly::ggplotly(p)

7.3.2.4 Code e8231bdf

load("data/7486fb94.RData")
library(ggplot2)
p=ggplot(d,aes(w1,t1))
p=p+xlab("Wind direction 7.30 am")
p=p+ylab("Temperature 7.30 am")
p=p+geom_jitter()
plotly::ggplotly(p)

load("data/7486fb94.RData")
library(ggplot2)
p=ggplot()
p=p+xlab("Wind direction 7.30 am")
p=p+ylab("Temperature 7.30 am")
p=p+geom_jitter(data=d,aes(w1,t1),size=2,color="blue")
plotly::ggplotly(p)

7.3.2.5 Code b4ee328a

load("data/7486fb94.RData")
library(ggplot2)
p=ggplot(d,aes(t1,t3))
p=p+geom_point(color="red",size=3)
p=p+xlab("Temperature 7.30 am")
p=p+ylab("Temperature 9.30 pm")
p=p+scale_x_continuous(breaks=seq(-5,25,5))
p=p+scale_y_continuous(breaks=seq(-5,25,5))
p=p+theme(axis.text=element_text(size=16),
          axis.title=element_text(size=16,face="bold"))
plotly::ggplotly(p)

7.3.2.6 Code 8076ef9d

load("data/7486fb94.RData")
library(ggplot2)
p=ggplot(d,aes(t1,t3))
p=p+geom_line(color="red",size=1)
p=p+xlab("Temperature 7.30 am")
p=p+ylab("Temperature 9.30 pm")
p=p+scale_x_continuous(breaks=seq(-5,25,5))
p=p+scale_y_continuous(breaks=seq(-5,25,5))
p=p+theme(axis.text=element_text(size=16),
          axis.title=element_text(size=16,face="bold"))
plotly::ggplotly(p)

7.3.2.7 Code fd54469e

load("data/7486fb94.RData")
library(ggplot2)
p=ggplot(d,aes(da1,t2))
p=p+geom_line(color="red",size=1)
p=p+xlab("Date (1990)")
p=p+ylab("Temperature 2.30 pm")
p=p+scale_x_date(date_labels = "%m-%d"
                 ,date_breaks = "1 week"
                 ,limit=c(as.Date("1990-01-01"),as.Date("1990-12-31"))
                 )
p=p+scale_y_continuous(breaks=seq(-5,35,5))
p=p+theme(axis.text=element_text(size=11),
          ,axis.title=element_text(size=16,face="bold")
          ,axis.text.x = element_text(angle=75)
          )
plotly::ggplotly(p)

7.3.2.8 Code c025a413

load("data/7486fb94.RData")
m=lm(t3~t1,d)
coef(m)
library(ggplot2)
p=ggplot(d,aes(t1,t3))
p=p+geom_point(color="red",size=3)
p=p+stat_function(fun=function(x) coef(m)[2]*x+coef(m)[1],color="black",size=3)
p=p+xlab("Temperature 7.30 am")
p=p+ylab("Temperature 9.30 pm")
p=p+scale_x_continuous(breaks=seq(-5,25,5))
p=p+scale_y_continuous(breaks=seq(-5,25,5))
p=p+theme(axis.text=element_text(size=16),
          axis.title=element_text(size=16,face="bold"))
plotly::ggplotly(p)

load("data/7486fb94.RData")
m=lm(t3~t1,d)
coef(m)

(Intercept)          t1 
  2.0786062   0.9625644

7.3.2.9 Code 22e19e4b TODO

predicted measured plot sd in abhängigkeit von t1,rfm usw.

7.3.2.10 Code c1885779

load("data/7486fb94.RData")
library(ggplot2)
# Colors:
#######################
library(pals)
cols1=as.character(polychrome())
ind=1:length(cols1)
cols1=c(cols1[ind[!ind %in% c(1,2,22)]],cols1[ind[!ind %in% c(1,2,22)]])

cols2=c("black","red","blue","green",  "yellow", "brown"    ,"magenta","cyan","darkgreen","orange","darkgray","hotpink","cornflowerblue","lawngreen","gold","chocolate1","red4","deeppink","navy","black")

cols=c(cols2,cols1)



p=ggplot()



leg=NULL
far=NULL
lab=NULL
i=0

### USER INPUT: specify x and y axis columns and column labels colyl
i=i+1
colx="da1"
coly="t1"
colyl="7.30 am"

# leg=append(leg,'colyl')
cmd=paste("leg=append(leg,'",colyl,"')",sep=""  )
eval(parse(text=cmd))

# far=append(far,c('colyl'='cols[i]')) 
cmd=paste("far=append(far,c('",colyl,"'='",cols[i],"'))",sep="")
eval(parse(text=cmd))

# lab=append(lab,c('colyl'='colyl'))
cmd=paste("lab=append(lab,c('",colyl,"'='",colyl,"'))",sep="")
eval(parse(text=cmd))

#  p=p+geom_line(data=d,aes(colx,coly,color='colyl'),size=1)
cmd=paste("p=p+geom_line(data=d,aes(",colx,",",coly,",color='",colyl,"'),size=1)",sep="")
eval(parse(text=cmd))



### USER INPUT: specify x and y axis columns and column labels colyl
i=i+1
colx="da1"
coly="t2"
colyl="2.30 pm"



# leg=append(leg,'colyl')
cmd=paste("leg=append(leg,'",colyl,"')",sep=""  )
eval(parse(text=cmd))

# far=append(far,c('colyl'='cols[i]')) 
cmd=paste("far=append(far,c('",colyl,"'='",cols[i],"'))",sep="")
eval(parse(text=cmd))

# lab=append(lab,c('colyl'='colyl'))
cmd=paste("lab=append(lab,c('",colyl,"'='",colyl,"'))",sep="")
eval(parse(text=cmd))

#  p=p+geom_line(data=d,aes(colx,coly,color='colyl'),size=1)
cmd=paste("p=p+geom_line(data=d,aes(",colx,",",coly,",color='",colyl,"'),size=1)",sep="")
eval(parse(text=cmd))




### USER INPUT: specify x and y axis columns and column labels colyl
i=i+1
colx="da1"
coly="t3"
colyl="9.30 pm"



# leg=append(leg,'colyl')
cmd=paste("leg=append(leg,'",colyl,"')",sep=""  )
eval(parse(text=cmd))

# far=append(far,c('colyl'='cols[i]')) 
cmd=paste("far=append(far,c('",colyl,"'='",cols[i],"'))",sep="")
eval(parse(text=cmd))

# lab=append(lab,c('colyl'='colyl'))
cmd=paste("lab=append(lab,c('",colyl,"'='",colyl,"'))",sep="")
eval(parse(text=cmd))

#  p=p+geom_line(data=d,aes(colx,coly,color='colyl'),size=1)
cmd=paste("p=p+geom_line(data=d,aes(",colx,",",coly,",color='",colyl,"'),size=1)",sep="")
eval(parse(text=cmd))



# Legend:
#######################
p=p+scale_color_manual(
    name = ""
   ,values=far
   ,labels=lab
   ,limits=leg
                       )
p=p+scale_x_date(date_labels = "%m-%d"
                 ,date_breaks = "1 week"
                 ,limit=c(as.Date("1990-01-01"),as.Date("1990-12-31"))
                 )


p=p+theme( axis.text.x = element_text(angle=75))


p=p+xlab("Date (1990)")
p=p+ylab("Temperature [°C]")


library(plotly)
# ggplotly(p)
ggplotly(p) %>%
layout(legend = list(orientation = "h", x = 0, y =1.1))

7.3.2.11 Code f5d8e726

# Show commands?
shcmd=T


load("data/7486fb94.RData")
library(ggplot2)
# Colors:
#######################
library(pals)
cols1=as.character(polychrome())
ind=1:length(cols1)
cols1=c(cols1[ind[!ind %in% c(1,2,22)]],cols1[ind[!ind %in% c(1,2,22)]])

cols2=c("black","red","blue","green",  "yellow", "brown"    ,"magenta","cyan","darkgreen","orange","darkgray","hotpink","cornflowerblue","lawngreen","gold","chocolate1","red4","deeppink","navy","black")

cols=c(cols2,cols1)



p=ggplot()



leg=NULL
far=NULL
lab=NULL
i=0

### USER INPUT: specify x and y axis columns and column labels colyl
colx="da1"
coly=c("t1","t2","t3")
colyl=c("7.30 am","2.30 pm","9.30 pm")

for (colya in coly) {
i=i+1

if (shcmd) cat("i:",i,"\n")

# leg=append(leg,'colyl')
cmd=paste("leg=append(leg,'",colyl[i],"')",sep=""  )
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

# far=append(far,c('colyl'='cols[i]')) 
cmd=paste("far=append(far,c('",colyl[i],"'='",cols[i],"'))",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

# lab=append(lab,c('colyl'='colyl'))
cmd=paste("lab=append(lab,c('",colyl[i],"'='",colyl[i],"'))",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

#  p=p+geom_line(data=d,aes(colx,coly,color='colyl'),size=1)
cmd=paste("p=p+geom_line(data=d,aes(",colx,",",colya,",color='",colyl[i],"'),size=1)",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

}


if (shcmd) far

if (shcmd) lab

if (shcmd) leg

# Legend:
#######################
p=p+scale_color_manual(
    name = ""
   ,values=far
   ,labels=lab
   ,limits=leg
                       )
p=p+scale_x_date(date_labels = "%m-%d"
                 ,date_breaks = "1 week"
                 ,limit=c(as.Date("1990-01-01"),as.Date("1990-12-31"))
                 )


p=p+theme( axis.text.x = element_text(angle=75))


p=p+xlab("Date (1990)")
p=p+ylab("Temperature [°C]")


library(plotly)
# ggplotly(p)
ggplotly(p) %>%
layout(legend = list(orientation = "h", x = 0, y =1.1))

# Show commands?
shcmd=T


load("data/7486fb94.RData")
library(ggplot2)
# Colors:
#######################
library(pals)
cols1=as.character(polychrome())
ind=1:length(cols1)
cols1=c(cols1[ind[!ind %in% c(1,2,22)]],cols1[ind[!ind %in% c(1,2,22)]])

cols2=c("black","red","blue","green",  "yellow", "brown"    ,"magenta","cyan","darkgreen","orange","darkgray","hotpink","cornflowerblue","lawngreen","gold","chocolate1","red4","deeppink","navy","black")

cols=c(cols2,cols1)



p=ggplot()



leg=NULL
far=NULL
lab=NULL
i=0

### USER INPUT: specify x and y axis columns and column labels colyl
colx="da1"
coly=c("t1","t2","t3")
colyl=c("7.30 am","2.30 pm","9.30 pm")

for (colya in coly) {
i=i+1

if (shcmd) cat("i:",i,"\n")

# leg=append(leg,'colyl')
cmd=paste("leg=append(leg,'",colyl[i],"')",sep=""  )
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

# far=append(far,c('colyl'='cols[i]')) 
cmd=paste("far=append(far,c('",colyl[i],"'='",cols[i],"'))",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

# lab=append(lab,c('colyl'='colyl'))
cmd=paste("lab=append(lab,c('",colyl[i],"'='",colyl[i],"'))",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

#  p=p+geom_line(data=d,aes(colx,coly,color='colyl'),size=1)
cmd=paste("p=p+geom_line(data=d,aes(",colx,",",colya,",color='",colyl[i],"'),size=1)",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

}

i: 1 
leg=append(leg,'7.30 am') 
far=append(far,c('7.30 am'='black')) 
lab=append(lab,c('7.30 am'='7.30 am')) 
p=p+geom_line(data=d,aes(da1,t1,color='7.30 am'),size=1) 
i: 2 
leg=append(leg,'2.30 pm') 
far=append(far,c('2.30 pm'='red')) 
lab=append(lab,c('2.30 pm'='2.30 pm')) 
p=p+geom_line(data=d,aes(da1,t2,color='2.30 pm'),size=1) 
i: 3 
leg=append(leg,'9.30 pm') 
far=append(far,c('9.30 pm'='blue')) 
lab=append(lab,c('9.30 pm'='9.30 pm')) 
p=p+geom_line(data=d,aes(da1,t3,color='9.30 pm'),size=1)

if (shcmd) far

7.30 am 2.30 pm 9.30 pm 
"black"   "red"  "blue"

if (shcmd) lab

  7.30 am   2.30 pm   9.30 pm 
"7.30 am" "2.30 pm" "9.30 pm"

if (shcmd) leg

[1] "7.30 am" "2.30 pm" "9.30 pm"

# Legend:
#######################
p=p+scale_color_manual(
    name = ""
   ,values=far
   ,labels=lab
   ,limits=leg
                       )
p=p+scale_x_date(date_labels = "%m-%d"
                 ,date_breaks = "1 week"
                 ,limit=c(as.Date("1990-01-01"),as.Date("1990-12-31"))
                 )


p=p+theme( axis.text.x = element_text(angle=75))


p=p+xlab("Date (1990)")
p=p+ylab("Temperature [°C]")


library(plotly)
# ggplotly(p)
ggplotly(p) %>%
layout(legend = list(orientation = "h", x = 0, y =1.1))

7.3.2.12 Code e7ffca3e

TODO Farben aus: mails/dominik-2023-05-24/

# Show commands?
shcmd=T # note: switched off after one iteration of the for loop

load("data/c55d0a7b.RData")

library(ggplot2)
# Colors:
#######################
library(pals)
cols1=as.character(polychrome())
ind=1:length(cols1)
cols1=c(cols1[ind[!ind %in% c(1,2,22)]],cols1[ind[!ind %in% c(1,2,22)]])

cols2=c("black","red","blue","green",  "yellow", "brown"    ,"magenta","cyan","darkgreen","orange","darkgray","hotpink","cornflowerblue","lawngreen","gold","chocolate1","red4","deeppink","navy","black")

cols=c(cols2,cols1)






leg=NULL
far=NULL
lab=NULL
i=0

### USER INPUT: specify x and y axis columns and column labels colyl
years=format(d$da1, format="%Y")
yearsl=levels(as.factor(years))
d$yday=as.numeric(format(d$da1, format="%j"))
colx="yday"
coly="t2"
colyl=yearsl



p=ggplot()


for (year in yearsl) {
i=i+1

dt=d[years==year,]

if (shcmd) cat("i:",i,"\n")

# leg=append(leg,'colyl')
cmd=paste("leg=append(leg,'",colyl[i],"')",sep=""  )
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

# far=append(far,c('colyl'='cols[i]')) 
cmd=paste("far=append(far,c('",colyl[i],"'='",cols[i],"'))",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

# lab=append(lab,c('colyl'='colyl'))
cmd=paste("lab=append(lab,c('",colyl[i],"'='",colyl[i],"'))",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

#  p=p+geom_line(data=d,aes(colx,coly,color='colyl'),size=1)
cmd=paste("p=p+geom_line(data=dt,aes(",colx,",",coly,",color='",colyl[i],"'),size=1)",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

shcmd=F

}


if (shcmd) far

if (shcmd) lab

if (shcmd) leg

# Legend:
#######################
p=p+scale_color_manual(
    name = ""
   ,values=far
   ,labels=lab
   ,limits=leg
                       )

p=p+theme( axis.text.x = element_text(angle=75))

p=p+scale_x_continuous(breaks=seq(0,360,10))
p=p+scale_y_continuous(breaks=seq(-20,40,5),limits = c(-15,45))

p=p+xlab("Day")
p=p+ylab("Temperature [°C]")


library(plotly)
# ggplotly(p)
ggplotly(p) %>%
layout(legend = list(orientation = "h", x = 0, y =1.1))

i: 1 
leg=append(leg,'1961') 
far=append(far,c('1961'='black')) 
lab=append(lab,c('1961'='1961')) 
p=p+geom_line(data=dt,aes(yday,t2,color='1961'),size=1)

7.3.2.13 Code dd4d414c

# Show commands?
shcmd=T # note: switched off after one iteration of the for loop

load("data/c55d0a7b.RData")

library(ggplot2)
# Colors:
#######################
library(pals)
cols1=as.character(polychrome())
ind=1:length(cols1)
cols1=c(cols1[ind[!ind %in% c(1,2,22)]],cols1[ind[!ind %in% c(1,2,22)]])

cols2=c("black","red","blue","green",  "yellow", "brown"    ,"magenta","cyan","darkgreen","orange","darkgray","hotpink","cornflowerblue","lawngreen","gold","chocolate1","red4","deeppink","navy","black")

cols=c(cols2,cols1)






leg=NULL
far=NULL
lab=NULL
i=0

### USER INPUT: specify x and y axis columns and column labels colyl
years=format(d$da1, format="%Y")
yearsl=levels(as.factor(years))
# d$yday=as.numeric(format(d$da1, format="%j"))
colx="w2"
coly="t2"
colyl=yearsl



p=ggplot()


for (year in yearsl) {
i=i+1

dt=d[years==year,]

if (shcmd) cat("i:",i,"\n")

# leg=append(leg,'colyl')
cmd=paste("leg=append(leg,'",colyl[i],"')",sep=""  )
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

# far=append(far,c('colyl'='cols[i]')) 
cmd=paste("far=append(far,c('",colyl[i],"'='",cols[i],"'))",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

# lab=append(lab,c('colyl'='colyl'))
cmd=paste("lab=append(lab,c('",colyl[i],"'='",colyl[i],"'))",sep="")
eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

#  p=p+geom_line(data=d,aes(colx,coly,color='colyl'),size=1)
# cmd=paste("p=p+geom_line(data=dt,aes(",colx,",",coly,",color='",colyl[i],"'),size=1)",sep="")
cmd=paste("p=p+geom_jitter(data=dt,aes(",colx,",",coly,",color='",colyl[i],"'),size=1)",sep="")

eval(parse(text=cmd))
if (shcmd) cat(cmd,"\n")

shcmd=F

}


if (shcmd) far

if (shcmd) lab

if (shcmd) leg

# Legend:
#######################
p=p+scale_color_manual(
    name = ""
   ,values=far
   ,labels=lab
   ,limits=leg
                       )

# p=p+theme( axis.text.x = element_text(angle=75))

# p=p+scale_x_continuous(breaks=seq(0,360,10))
p=p+scale_y_continuous(breaks=seq(-20,40,5),limits = c(-15,45))

p=p+xlab("Wind direction 2.30 pm")
p=p+ylab("Temperature 2.30 pm [°C]")


library(plotly)
# ggplotly(p)
ggplotly(p) %>%
layout(legend = list(orientation = "h", x = 0, y =1.0))

i: 1 
leg=append(leg,'1961') 
far=append(far,c('1961'='black')) 
lab=append(lab,c('1961'='1961')) 
p=p+geom_jitter(data=dt,aes(w2,t2,color='1961'),size=1)

7.3.2.14 Code 542ad9cc TODO

3D-Grafiken aus crashr

Plot einer 3D-Funktion

7.3.3 Programming

7.3.3.1 Functions

7.3.3.1.1 Code 0cc6a21b

f=function(x) x^2
f(2)
curve(f,-1,1)
f=function(x,y) x^2+y^2
f(2,3)

f=function(x) x^2
f(2)

[1] 4

curve(f,-1,1)

f=function(x,y) x^2+y^2
f(2,3)

[1] 13

7.3.3.2 Repeated operations

7.3.3.2.1 Code 775f4042

i=28
cat("Der Wert von i ist:",i)
for (i in 1:3) cat("Der Wert von i ist:",i,"\n")
for (i in 1:3) cat("Der Wert von i ist:",i)
sum=0
for (i in 1:3) sum=sum+i
sum
sum=0
for (i in 1:3) {
  cat("Der Wert von i ist:",i,"\n")
  cat("Der Wert von sum ist:",sum,"\n")
  sum=sum+i
  cat("Der Wert von sum=sum+i ist:",sum,"\n")
  cat("\n")
 }
for (i in c(1,3,5)) cat("Der Wert von i ist:",i,"\n")
for (col in c("t1","t2","t3")) cat("Der Wert von col ist:",col,"\n")
load("data/7486fb94.RData")
for (col in c("t1","t2","t3")){
  mw=mean(d[,col])
  cat("Der Mittelwert der Spalte",col,"ist:",mw,"\n")
}
load("data/7486fb94.RData")
for (col in colnames(d)){
  cat("Datentyp der Spalte",col,"ist:",class(d[,col]),"\n")
}
sapply(1:3,function(x) x^2)
load("data/7486fb94.RData")
sapply(d,class)

i=28
cat("Der Wert von i ist:",i)

Der Wert von i ist: 28

for (i in 1:3) cat("Der Wert von i ist:",i,"\n")

Der Wert von i ist: 1 
Der Wert von i ist: 2 
Der Wert von i ist: 3

for (i in 1:3) cat("Der Wert von i ist:",i)

Der Wert von i ist: 1Der Wert von i ist: 2Der Wert von i ist: 3

sum=0
for (i in 1:3) sum=sum+i
sum

[1] 6

sum=0
for (i in 1:3) {
  cat("Der Wert von i ist:",i,"\n")
  cat("Der Wert von sum ist:",sum,"\n")
  sum=sum+i
  cat("Der Wert von sum=sum+i ist:",sum,"\n")
  cat("\n")
 }

Der Wert von i ist: 1 
Der Wert von sum ist: 0 
Der Wert von sum=sum+i ist: 1 

Der Wert von i ist: 2 
Der Wert von sum ist: 1 
Der Wert von sum=sum+i ist: 3 

Der Wert von i ist: 3 
Der Wert von sum ist: 3 
Der Wert von sum=sum+i ist: 6

for (i in c(1,3,5)) cat("Der Wert von i ist:",i,"\n")

Der Wert von i ist: 1 
Der Wert von i ist: 3 
Der Wert von i ist: 5

for (col in c("t1","t2","t3")) cat("Der Wert von col ist:",col,"\n")

Der Wert von col ist: t1 
Der Wert von col ist: t2 
Der Wert von col ist: t3

load("data/7486fb94.RData")
for (col in c("t1","t2","t3")){
  mw=mean(d[,col])
  cat("Der Mittelwert der Spalte",col,"ist:",mw,"\n")
}

Der Mittelwert der Spalte t1 ist: 8.723014 
Der Mittelwert der Spalte t2 ist: 14.30932 
Der Mittelwert der Spalte t3 ist: 10.47507

load("data/7486fb94.RData")
for (col in colnames(d)){
  cat("Datentyp der Spalte",col,"ist:",class(d[,col]),"\n")
}

Datentyp der Spalte da ist: character 
Datentyp der Spalte p1 ist: numeric 
Datentyp der Spalte p2 ist: numeric 
Datentyp der Spalte p3 ist: numeric 
Datentyp der Spalte pd ist: numeric 
Datentyp der Spalte tmax ist: numeric 
Datentyp der Spalte tmin ist: numeric 
Datentyp der Spalte emin ist: numeric 
Datentyp der Spalte t1 ist: numeric 
Datentyp der Spalte t2 ist: numeric 
Datentyp der Spalte t3 ist: numeric 
Datentyp der Spalte rm1 ist: integer 
Datentyp der Spalte rm2 ist: integer 
Datentyp der Spalte rm3 ist: integer 
Datentyp der Spalte rmd ist: integer 
Datentyp der Spalte w1 ist: character 
Datentyp der Spalte w2 ist: character 
Datentyp der Spalte w3 ist: character 
Datentyp der Spalte ws1 ist: integer 
Datentyp der Spalte ws2 ist: integer 
Datentyp der Spalte ws3 ist: integer 
Datentyp der Spalte wsd ist: numeric 
Datentyp der Spalte cc1 ist: integer 
Datentyp der Spalte cc2 ist: integer 
Datentyp der Spalte cc3 ist: integer 
Datentyp der Spalte ccd ist: numeric 
Datentyp der Spalte sd ist: numeric 
Datentyp der Spalte vr1 ist: character 
Datentyp der Spalte vr2 ist: character 
Datentyp der Spalte vr3 ist: character 
Datentyp der Spalte pr1 ist: numeric 
Datentyp der Spalte prt1 ist: character 
Datentyp der Spalte pr2 ist: numeric 
Datentyp der Spalte prt2 ist: character 
Datentyp der Spalte pr3 ist: numeric 
Datentyp der Spalte prt3 ist: character 
Datentyp der Spalte pd.1 ist: numeric 
Datentyp der Spalte pdt ist: character 
Datentyp der Spalte sh ist: integer 
Datentyp der Spalte wsdmax ist: numeric 
Datentyp der Spalte da1 ist: Date

sapply(1:3,function(x) x^2)

[1] 1 4 9

load("data/7486fb94.RData")
sapply(d,class)

         da          p1          p2          p3          pd        tmax 
"character"   "numeric"   "numeric"   "numeric"   "numeric"   "numeric" 
       tmin        emin          t1          t2          t3         rm1 
  "numeric"   "numeric"   "numeric"   "numeric"   "numeric"   "integer" 
        rm2         rm3         rmd          w1          w2          w3 
  "integer"   "integer"   "integer" "character" "character" "character" 
        ws1         ws2         ws3         wsd         cc1         cc2 
  "integer"   "integer"   "integer"   "numeric"   "integer"   "integer" 
        cc3         ccd          sd         vr1         vr2         vr3 
  "integer"   "numeric"   "numeric" "character" "character" "character" 
        pr1        prt1         pr2        prt2         pr3        prt3 
  "numeric" "character"   "numeric" "character"   "numeric" "character" 
       pd.1         pdt          sh      wsdmax         da1 
  "numeric" "character"   "integer"   "numeric"      "Date"

7.3.3.3 Conditional operations

7.3.3.3.1 Code 001863d3

if (FALSE) cat("Hello world! \n")
if (TRUE) cat("Hello world! \n")
if (TRUE) 
  cat("The value is TRUE. \n")  else 
  cat("The value is FALSE. \n")
if (FALSE) 
  cat("The value is TRUE. \n")  else 
  cat("The value is FALSE. \n")
a=-1
if (a<0) 
  cat("a is negative. \n")  else 
  cat("a is positive. \n")
a=1
if (a<0) 
  cat("a is negative. \n")  else 
  cat("a is positive. \n")
a=1;b=1
if (a>0 & b>0) 
  cat("a,b positive. \n")  else 
  if (a<0 & b<0) cat("a,b negative. \n") else
  cat("Something else. \n") 
a=-1;b=-1
if (a>0 & b>0) 
  cat("a,b positive. \n")  else 
  if (a<0 & b<0) cat("a,b negative. \n") else
  cat("Something else. \n") 
a=-1;b=0
if (a>0 & b>0) 
  cat("a,b positive. \n")  else 
  if (a<0 & b<0) cat("a,b negative. \n") else
  cat("Something else :-) \n")

if (FALSE) cat("Hello world! \n")
if (TRUE) cat("Hello world! \n")

Hello world!

if (TRUE) 
  cat("The value is TRUE. \n")  else 
  cat("The value is FALSE. \n")

The value is TRUE.

if (FALSE) 
  cat("The value is TRUE. \n")  else 
  cat("The value is FALSE. \n")

The value is FALSE.

a=-1
if (a<0) 
  cat("a is negative. \n")  else 
  cat("a is positive. \n")

a is negative.

a=1
if (a<0) 
  cat("a is negative. \n")  else 
  cat("a is positive. \n")

a is positive.

a=1;b=1
if (a>0 & b>0) 
  cat("a,b positive. \n")  else 
  if (a<0 & b<0) cat("a,b negative. \n") else
  cat("Something else. \n")

a,b positive.

a=-1;b=-1
if (a>0 & b>0) 
  cat("a,b positive. \n")  else 
  if (a<0 & b<0) cat("a,b negative. \n") else
  cat("Something else. \n")

a,b negative.

a=-1;b=0
if (a>0 & b>0) 
  cat("a,b positive. \n")  else 
  if (a<0 & b<0) cat("a,b negative. \n") else
  cat("Something else :-) \n")

Something else :-)