---
title: "R Notebook"
output: html_notebook
---


```{r echo=FALSE, message=FALSE, warning=FALSE}


install.packages(pkgs = "")
install.packages(file_name_and_path, repos = NULL, type="source")


# RStudio - cran

# RStudio - archiv

# GitHub
library(devtools)
install_github("Momocs",username="jfpalomeque")
install_github("jfpalomeque/Momocs")


# Bioconductor
# https://www.bioconductor.org/
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.14")

BiocManager::available()

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("Rchemcpp")
BiocManager::install("biodb")


# Search

library(pkgsearch)
pkg_search_addin(query = "", viewer = c("dialog", "browser"))

# https://www.itl.nist.gov/div898/education/datasets.htm


library(NCmisc)
require(NCmisc)

summarise.r.datasets()
summarise.r.datasets(filter=TRUE,"matrix")

```


Přiřazení hodnoty proměnné

```{r echo=FALSE, warning=FALSE}

x <- -7
x

x = -7; x

(x = -7)


library(schoolmath)
schoolmath::is.negative(x)
schoolmath::is.positive(x)
schoolmath::is.real.positive(x)

abs(x)

```


Sčítání, odčítání, násobení a dělení

```{r echo=FALSE, message=FALSE, warning=FALSE}

3+3
9-3
3*3
9/3

x = 9
x+3
x-3
x*3
x/3
(x-3)/(x/3)

1/0
Inf*Inf
Inf/Inf
0/0


# celociselne deleni (intager division)
13%/%2
library(numbers)
div(13,2)  

# zbytek po deleni (reminder after division)
13%%2
library(numbers)
mod(13,2) 

```


Převod desetinných čísel na zlomky

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(MASS)
fractions(-0.1666667)
fractions(0.14)
fractions(0.4)

library(schoolmath)
decimal2fraction(0.1,6)

library(fractional)
fractional(-0.1666667, eps = 1e-06, maxConv = 20, sync = FALSE)
numerators(-0.1666667) # vypise citatele (numerator)
denominators(-0.1666667) # vypise jmenovatele (denominator)
numerical(-1/6)
unfractional(-1/6)

# zjednodusení zlomku (fractions simplification)
library(schoolmath)
cancel.fraction(42, 56)
cancel.fraction(-167, 100)


# Nejvetsi spolecny delitel (greatest common divisor)
x = 12
y = 8
library(schoolmath)
schoolmath::gcd(x, y)
library(DescTools)
DescTools::GCD(x, y, na.rm = FALSE)
DescTools::GCD(c(12,10,8), na.rm = FALSE)


# Nejmensi spolecny nasobek (smallest common multiple)
x = 12
y = 8
library(schoolmath)
schoolmath::scm(x, y)
library(DescTools)
DescTools::LCM(x, y)
DescTools::LCM(c(12,10,8), na.rm = FALSE)

```


Mocniny a odmocniny

```{r echo=FALSE, message=FALSE, warning=FALSE}

x = 9
y=4

sqrt(x)
x^(1/2)
x^(-1/2)

x**2
x^2
x^(-2)


# zda je cislo mocninou celeho cisla

library(numbers)
isIntpower(144)
isIntpower(8)

```


Exponenty a logaritmy

```{r echo=FALSE, message=FALSE, warning=FALSE}

### Exponenty ###
x = 2
10^x # exponent
10**x # exponent

x = 2
exp(x)   # e^x


### Prirozeny logaritmus ###
x = 10
log(x)  # base e
log(x, base = exp(1))  # base e
log(0)


### Dekadicky logaritmus ###
x = 10
log10(x) # base 10
log(x,base = 10)
log10(0)

```


Goniometrické funkce

```{r echo=FALSE, message=FALSE, warning=FALSE}

x = pi/2
cos(x)
sin(x)
tan(x)

x = 1
acos(x)
asin(x)
atan(x)


# prevod uhly vs radiany
library(pracma)
deg2rad(120) 
rad2deg(3.14)

```


Vypočítejte pH 0.001 M kyseliny chlorovodíkové, 0.01 M hydroxidu sodného, 0.15 M kyseliny octové (Ka = 1.75e-5) a 0.01 M amoniaku (Kb = 1.8e-5)

```{r echo=FALSE, message=FALSE, warning=FALSE}

# HCl
pH = -log10(0.001); pH


# NaOH
pH = 14 - (-log10(0.01)); pH


# CH3COOOH
pH = (-log10(1.75e-5)-log10(0.15))/2; pH

library(ch)
monoa(ka=1.75e-5, c=0.15, digits = 2)
mono(ka=1.75e-5, c=0.15, digits = 4, acid = TRUE, kw = 1e-14)


# NH4OH
pH = 14 - (-log10(1.8e-5)-log10(0.01))/2; pH

library(ch)
monob(ka=1.8e-5, c=0.01, digits = 2)
mono(ka=1.8e-5, c=0.01, digits = 4, acid = FALSE, kw = 1e-14)

```


Matematické konstanty

```{r echo=FALSE, message=FALSE, warning=FALSE}

### Ludolfovo cislo ####
pi
print(pi, digits=20)

### Eulerovo cislo ####
exp(1)
print(exp(1), digits=20) 

```


Vektory - generování řady čísel

```{r echo=FALSE, message=FALSE, warning=FALSE}

1:6
c(2:12)
rev(c(2:12))

seq(from=1, to=20, by=2)
seq(5, -5, -1)
seq(from=5, by=-1, along.with = 1:20)

rep(1, len = 9)
rep(1:4, len = 9)
rep(1:4, 2)
rep(1:4, each = 2)       # not the same.
rep(1:4, each = 2, len = 4)    # first 4 only.
rep(1:4, each = 2, times = 3)  

```


Vektory - řada náhodných čísel z daného rozmezí

```{r echo=FALSE, message=FALSE, warning=FALSE}

# set.seed(123)
sample(12:99, 9, replace = FALSE) # bez opakovani
sort(sample(1:20, 9, replace = FALSE))
sample(12:99, 9, replace = TRUE) # s opakovanim
sort(sample(1:20, 9, replace = TRUE))

rnorm(n=6, mean=0, sd=1)
rnorm(6,9,1.5)

rbinom(30,c(0,1),0.5) #mince
rbinom(30,c(1:6),(1/6)) #kostka (dice)

runif(n=5, min=0, max=1) # rovnomerne rozdeleni
runif(6,-2,2)

```


Indexy

```{r echo=FALSE, message=FALSE, warning=FALSE}
set.seed(123) # generator nahodnych cisel
x = sample(12:99,19, replace = FALSE)
is.vector(x)
x
x[1]
x[1:5]
x[c(1:5)]
x[length(x)]
x[(length(x)-1)]
x[-1]
x[-(1:5)]
x[-c(1:5)]
x[-length(x)]
x[(length(x)-3):length(x)]#posledni 4 hodnoty
x[c(1,4,7)]#vybrane hodnoty
which.max(x)
which.min(x)
library(Rmpfr)
which.min(x)
which(x<=50)
x[which(x<=50)]
x[x<=50]
which(x>=50)
x[which(x>=50)]
which(x == max(x))#maximalni hodnota
which(x == 10)
which(x != 10)
x[x > 40]#hodnoty vyssi nez 40.
x[x < 40]#hodnoty mensi nez 40
which(x > 10)#hodnoty vyssi nez 10.
which(x < 10)#hodnoty mensi nez 10.

```


Uspořádání vektoru

```{r echo=FALSE, message=FALSE, warning=FALSE}

x = sample(12:99,19, replace = FALSE)

rank(x)  # vhodne pro vektory, nevhodne pro matice 

sort(x, decreasing = FALSE)
sort(x, decreasing = TRUE)
sort(x[which(x>=50)])
order(x)  # poradi (indexy) pro sort
x[order(x)]

```


Množiny

```{r echo=FALSE, message=FALSE, warning=FALSE}

x = 10

x == 0
x != 0
x == 10
x != 10
x>0
x<0
x>=0
x>=10
x<=10
-x>=10

library(extraoperators)
# https://joshuawiley.com/extraoperators/
x%l%10 # <
x%le%10 # <=
x%g%10 # >
x%ge%10 # >=
 
 
x <- c(sort(sample(1:20, 3)))

x>5 & x<10 # 5 < x AND x < 11
x>5 | x<10 # 5 < x OR x < 11

library(extraoperators)
# https://joshuawiley.com/extraoperators/

5%gl%10 # # 5 < x AND x < 11 (5, 11)
5%gel%10 # 5 <= x AND x < 11 [5, 11)
5%gle%10 # 5 < x AND x <= 11 (5, 11]
5%gele%10 # 5 <= x AND x <= 11 [x, y]


y <- c(sort(sample(3:23, 7)))

union(x, y) # sjednoceni
intersect(x, y) # prunik
setdiff(x, y) # rozdil mnozin
setdiff(y, x) # rozdil mnozin


# kladne vs. zaporne cislo
x <- c(-1, -2, 3.02, 4, -5.2, 6, -7, 0) 
schoolmath::is.negative(x)
schoolmath::is.positive(x)
schoolmath::is.real.positive(x)

# sude vs. liche cislo
x <- c(1, 2, 3, 4, 5, 6, 7) 
schoolmath::is.even(x)
schoolmath::is.odd(x)
x %% 2 == 0 
x %% 2 != 0

```


Vektorová aritmetika

```{r}
a = c(1, 3, 5, 7) 
b = c(1, 2, 4, 8) 
5 * a
a + b 
a - b 
a * b 
a / b

u = c(10, 20, 30) 
v = c(1, 2, 3, 4, 5, 6, 7, 8, 9) 
u + v 

diff(u)
diff(v)


max(x)
min(x)
range(x)
library(DescTools)
Range(x)
sum(x)
cumsum(x)
prod(x)
cumprod(x)
diff(x)
sign(x)

```


Filtrování hodnot ve vektoru

```{r}

x = sample(12:99, 9, replace = TRUE)
x
unique(x)
table(x)
summary(as.factor(x))
# frekvence vyskytu daneho cisla
sum(x==32)
length(which(x==57))


numbers =c(4,23,4,23,5,43,54,56,657,67,67,435,453,435,7,65,34,435)
table(numbers)[2]==1
numbers[2]
sum(numbers==435)

duplicated(numbers) # identifying duplicated elements 
numbers[duplicated(numbers)]
unique(numbers) # extracting unique elements,

```


Matice

```{r}

x = sample(12:99,18, replace = FALSE); x
matrix(x,nrow=6,ncol=3,byrow = TRUE)
matrix(x,nrow=6,ncol=3,byrow = FALSE)


x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)

x12 = rbind(x1,x2)
x12 = cbind(x1,x2)
is.matrix(x12)
nrow(x12)
ncol(x12)

x12[order(x1),]
x12[order(x2),]
x12[,1]
x12[,"x"]
x12$x

colnames(x12)
colnames(x12)=c("X","Y")
rownames(x12)
rownames(x12) = LETTERS[1:length(x1)] 

```


Dataframes

```{r echo=FALSE, message=FALSE, warning=FALSE}

x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
z = LETTERS[1:length(x1)] 

x12 = rbind(x1,x2,z)
x12 = cbind(x1,x2,z)
nrow(x12)
ncol(x12)
is.data.frame(x12)
is.matrix(x12)

x12 = as.data.frame(x12)
x12[,1] = as.numeric(x12[,1])
x12[,2] = as.numeric(x12[,2])

is.data.frame(x12)
colnames(x12)
rownames(x12)
as.matrix(x12)


mean(x1)
mean(as.numeric(x12[,1]))

```


Seznamy (lists)

```{r echo=FALSE, message=FALSE, warning=FALSE}

x= list(1,2,3)
names(x) = c("x1","x2","x3")
x$x1[1]
x[[2]][1]

x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
x3 = sample(12:99,19, replace = TRUE)
z = LETTERS[1:length(x1)] 

x123 = list(x1,x2,x3,z)
is.list(x123)
names(x123) = c("x1","x2","x3","z") 
x123$x1
x123[[1]]


x1 = sample(12:99,19, replace = FALSE)
x2 = sample(12:99,19, replace = TRUE)
x3 = sample(12:99,19, replace = TRUE)
LL = list(x1,x2,x3)
names(LL) = c("x1","x2","x3")

LL[[2]][3]
LL$x2[3]

```



Vypočítejte molekulové hmotnosti alkanů C1 - C12.

```{r echo=FALSE, message=FALSE, warning=FALSE}

HC = function(n){n*12.011 + (2*n + 2)*1.008}

hc112 = HC(c(1:12))
names(hc112) = c(1:12)
hc112

hc112 = sapply(c(1:12),HC)
names(hc112) = c(1:12)
hc112

hc112 = lapply(c(1:12),HC)
names(hc112) = c(1:12)
hc112

```


Vliv 3 ůzných druhů krmiva na přírůstek živé váhy (v kg) prasat, rozdělených do 4 různých skupin. 

```{r echo=FALSE, message=FALSE, warning=FALSE}

A = c(7.0, 16.0, 10.5, 13.5)
B = c(14.0, 15.5, 15.0, 21.0)
C = c(8.5, 16.5, 9.5, 13.5)
ps =cbind(A,B,C)
colnames(ps) = c("A","B","C")
rownames(ps) = c("I","II","III","IV")
ps

colSums(ps)
rowSums(ps)

colMeans(ps)
rowMeans(ps)

apply(ps,2,mean)
apply(ps,1,mean)

# deleni hodnot ve sloupci sumou sloupce
sup = apply(ps,2,sum)
sapply(c(1:3),function(i){ps[,i]/sup[i]})


library(reshape2)
melt(ps)
pr = setNames(melt(ps), c('serie','krmeni','prirustek'))
pr

as.data.frame(pr)
as.matrix(pr)

acast(pr,serie ~ krmeni)


st = stack(as.data.frame(ps)); st
unstack(st)


prirustek = pr[,"prirustek"]
serie = as.factor(pr[,"serie"])
krmeni = as.factor(pr[,"krmeni"])

aggregate(prirustek,list(krmeni), FUN=mean)
aggregate(prirustek,list(serie,krmeni), FUN=mean)

tapply(prirustek, list(krmeni), mean)
tapply(pr$prirustek, pr$krmeni, mean)



lps = list(A,B,C)
lps
names(lps) = c("A","B","C")
lps

unlist(lps)

do.call(cbind,lps)

sapply(lps,mean)
lapply(lps,mean)


st = stack(lps); st
unstack(st)

library(reshape2)
melt(lps)
lpr = setNames(melt(lps), c('krmeni','prirustek'))
lpr


```


Pipe operator (%>%) v R

```{r}

a = 3.14159
b = seq(from = a,10,3)
round(b,3)

a = 3.14159
seq(round(a,3),10,3)

# pipe operator
library(magrittr)
a = 3.14159
a %>% seq(10,3) %>% round(3)

```


Příkaz if-else

```{r}

x = sample(12:99,18, replace = FALSE); x
ifelse(x<50,0,1)
x[x<50] = 0
x[x>=50] = 1
x

```


Vypocítejte vazebnou energii atomového jádra pomoci Bethe - Weiszackerovy rovnice B = 14.0 - 13.1*A^(2/3) + 0.585*Z*(Z-1)/A^(1/3) - (18.1*(A-2*Z)^2)/A + C/A, kde
A - nukleonove cislo, Z - atomove cislo, C - konstanta, pro jádra se sudým počtem protonů i neutronů je rovna 132, 
pro jádra s lichým počtem protonů i neutronů je rovna -132, pro lichy počet protonů a sudý počet neutronů a naopak je rovna 0.

```{r}

# alpha particle
A = 4 
Z = 2

# U-235
A = 235 
Z = 92

# funkce
N = A-Z
if (N%%2==0 & Z%%2==0){C = 132} else if (N%%2!=0 & Z%%2!=0){C = -132} else {C = 0}
B = 14.0 - 13.1*A^(2/3) + 0.585*Z*(Z-1)/A^(1/3) - (18.1*(A-2*Z)^2)/A + C/A
B # [MeV]
B/A # energie na 1 nukleon [MeV]

# polomer jadra: R = R0 * A^(1/3)
R0 = 1.4e-13 # [cm]
R = R0 * A^(1/3)
R = R*10e-15 
R   # [m]

```


Cyklus typu while

Vypočítejte molekulové hmotnosti alkanů C1 - C12 (while loop)

```{r}

n = 1
Mhc = NULL
while(n <= 12){
  hc = n*12.011 + (2*n + 2)*1.008
  Mhc = c(Mhc,hc)
  n = n+1
}
names(Mhc) = c(1:12)
Mhc

```


Cyklus typu for (for loop)

Vypočítejte molekulové hmotnosti alkanů C1 - C12

```{r}

Mhc = NULL
for(n in c(1:12)){
  hc = n*12.011 + (2*n + 2)*1.008
  Mhc = c(Mhc,hc) 
}
names(Mhc) = c(1:12)
Mhc

```


Vypočítejte střední kvadratickou rychlost molekul daných plynů v rozmezí teplot od 150 K do 500 K s krokem 50 K.

```{r}

el = c("H2","He", "O2", "Kr")
M = c(2.0159, 4.0026, 31.9988, 83.8) # [g / mol]
Ts = seq(from=150, to=500, by=50) # [K]
R = 8314.34 # [J / kmol K]
MKV = NULL
for(tt in Ts){
  mkv= rep(0, length(M))
  for(ii in c(1:length(M))){
    mkv[ii] = (3*R*tt/M[ii])^(1/2)
  }
  MKV = rbind(MKV, mkv)   
}
colnames(MKV) = el
rownames(MKV) = Ts
MKV
plot(0,0,xlim=c(min(Ts),max(Ts)),ylim = c(min(MKV),max(MKV)), type ="n",xlab="T [K]",ylab="str. kv. rychlost [m/s]")
for(ii in c(1:length(M))){points(Ts,MKV[,ii], type="p",col=ii,pch=16)}
legend("topleft",legend=el,col=c(1:length(M)),pch=16,cex=.8)

```


Fyzikální konstanty

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(marelac)
data(Constants)
Constants

# Univerzalni plynova konstanta
Constants$gasCt1
Constants$gasCt2
Constants$gasCt3
R = as.numeric(Constants$gasCt2[1]); R # hodnota
Constants$gasCt2[2] # jednotka


library(constants)
codata
syms_with_errors
syms_with_units

lookup("planck", ignore.case=TRUE)
lookup("avogadro", ignore.case=TRUE)
lookup("faraday", ignore.case=TRUE)
lookup("molar", ignore.case=TRUE)
lookup("molar mass", ignore.case=TRUE)
lookup("light", ignore.case=TRUE)
lookup("boltzmann", ignore.case=TRUE)
lookup("mass", ignore.case=TRUE)
lookup("mass constant", ignore.case=TRUE)
lookup("proton mass", ignore.case=TRUE)
lookup("electron mass", ignore.case=TRUE)

# Univerzalni plynova konstanta
lookup("gas", ignore.case=TRUE)
codata[134,]
R = as.numeric(codata[134,5]) # hodnota
Ru = as.numeric(codata[134,6]) # nejistota
codata[134,7] # jednotka

```


```{r}

library(ch)
period_table()

library(DescTools)
data(d.periodic)
d.periodic

library(loon.data)
data(elements)
elements

library(PeriodicTable)
data(periodicTable)
periodicTable
#
el = "Se"
pel = periodicTable[which(periodicTable$symb==el),]
pel$config
pel$mass
pel$Eneg
vel = pel$group-10  # pocet valencnich elektronu atomu

library(ChemmineR) # bioc
data(atomprop)
atomprop

library(CHNOSZ)
thermo()$element
thermo()$element[,c(1,4)]


#### Atomove hmotnosti

library(marelac)
data(AtomicWeight)
AtomicWeight

AtomicWeight[AtomicWeight$Symbol == "C",]
atomicweight$C 

molweight("H") 
molweight("C")

library(CHNOSZ)
AW = CHNOSZ::mass("Ca"); AW

library(biogas)
molMass("C")
molMass('Fe')

```


Určete délku měděného drátu o průměru 1.6 mm jehož hmotnost je 0.8 kg.

```{r}

library(PeriodicTable)
data(periodicTable)
ptCu = periodicTable[which(periodicTable$symb=="Cu"),]
ptCu$density

hm = 0.8 # kg
rho = ptCu$density*1000 # kg/m3
# rho = 8900 # kg/m3
d = 1.6e-3 # m

V = hm/rho
S = pi*(d/2)^2
L = V/S; L

paste("Délka drátu je", round(L,0), "m.")

```


Molekulové hmotnosti

```{r}

library(enviPat)
data(chemforms)
chemforms
formula1<-c("C12H10Cl10")
formula2<-c("C3H5")
# odecitani 2 vzorcu
subform(formula1,formula2)
# scitani 2 vzorcu
mergeform(formula1,formula2)
# nasobeni vzorce
formula1<-c("C12H10Cl10")
multiform(formula1,3)


library(marelac)
molweight("CO2") 
molweight("HCO3") 
molweight("H3PO4")
molweight("CH3CH2CHCHCH2CHCHCH2CHCHCH2CHCHCH2CHCH(CH2)3COOH") # eicosapentaenoic acid (EPA) 
molweight("C20H30O2")
molweight(c("C2H5OH", "CO2", "H2O")) 
molweight(c("SiOFH4", "NaHCO3", "C6H12O6", "Ca(HCO3)2", "Pb(NO3)2", "(NH4)2SO4"))


library(CHNOSZ)
MW = CHNOSZ::mass("CaCO3"); MW
# pocet atomu prvku ve vzorci
ce = count.elements("CaCO3"); ce
as.chemical.formula(ce, drop.zero = TRUE)

me = makeup("C6H12O6", multiplier = 1, sum = TRUE, count.zero = TRUE); me
ce = count.elements("C6H12O6"); ce
i2A("C6H12O6")
as.chemical.formula(ce, drop.zero = TRUE)
as.chemical.formula(me, drop.zero = TRUE)


library(biogas)
molMass("C6H12O6")
molMass(c("C6H12O6", "CH3COOH"))
molMass("H3C(CH2)5COOH")
molMass('FeSO4(H2O)7') # hydrates
molMass("(C6H12O6)0.24999 (H3COOH)0.75001")


library(BRAIN) # Bioc
calculateAverageMass(list(C=60, H=30, N=0, O=30, S=0))


library(MSbox)
MSbox::mass('C7H7O')
MSbox::mass('C7h7o1')
MSbox::mass('C7H7O', caseSensitive = TRUE)

```


Uhlovodík o molekulové hmotnosti 84 g/mol obsahuje 85.71 % uhlíku a 14.29 % vodíku. Určete molekulový vzorec.

```{r}

M = 84 # g/mol
mc = 0.8571
mh = 0.1429
library(marelac)
mC = molweight("C") # g/mol
mH = molweight("H") # g/mol

rs = c(mc/mC,mh/mH)/min(mc/mC,mh/mH)
rs = round(rs,0)
n = M/sum(rs*c(mC,mH))
n = round(n,0)
form = rs*n; form
gsub(" ", "", paste("C",form[1],"H",form[2]))


library(Rdisop)
dmp = decomposeMass(mass=M, ppm=1, mzabs=0.1,  elements=initializeElements(c("C","H")), z=1, minElements="C4", maxElements="C10")
dmp$formula

library(rcdk)
generate.formula(M, window=5, validation=FALSE, charge=1,elements=list(c("C",2,10),c("H",4,30)))

```


Izotopy

```{r}

library(enviPat)
data(isotopes)
isotopes

library(ecipex)
nistiso
nistiso$mass
nistiso$abundance
nistiso$element

library(IsoSpecR)
data(isotopicData)
isotopicData

library(CIAAWconsensus)
ciaaw.mass.2003
ciaaw.mass.2012
ciaaw.mass.2016

```


Izotopové patterny

```{r echo=FALSE, message=FALSE, warning=FALSE}


library(OrgMassSpecR)
IsotopicDistribution(formula = list(C = 1, H = 1, Cl = 3), charge = 1)
# allowed elements: C, H,N, O, S, P, Br, Cl, F, Si.
# relative atomic masses of the elements are from the NIST Physical Reference Data 
# https://www.nist.gov/pml/atomic-weights-and-isotopic-compositions-relative-atomic-masses.

mx <- IsotopicDistribution(formula = list(C = 9, H = 4, Br = 5, Cl = 1, N = 2))
library(lattice)
xyplot(percent ~ mz, data = mx,type = "h",xlab = "m/z",ylab = "intensity (%)",main = "Isotopic Distribution, C9H4Br5ClN2")


# C60H30O30

library(OrgMassSpecR)
x <- IsotopicDistribution(formula = ListFormula("C60H30O30"))
plot(x$mz, x$percent,type = "h",xlab = "m/z",ylab = "intensity (%)",main = "Isotopic Distribution, C60H30O30",lwd=2)
abline(h=0,lty=2)

library(ecipex)
xx <- ecipex("C60H30O30", isoinfo = nistiso, limit = 0.0005,id = FALSE, sortby = "mass")
lmass = xx$C60H30O30$mass
labun = xx$C60H30O30$abundance
plot(lmass,labun,type="h",xlab="m/z",ylab="abundance")
abline(h=0,lty=2)

library(rcdk)
formula1 <- get.formula("C60H30O30", charge = 1)
isot1 <- get.isotopes.pattern(formula1,minAbund=0.00005)
plot(isot1[,1], isot1[,2],type = "h",xlab = "m/z",ylab = "intensity (%)",main = "Isotopic Distribution, C60H30O30",lwd=2)
abline(h=0,lty=2)

library(enviPat)
data(isotopes)

patt = isopattern(isotopes, "C60H30O30", threshold = 0.001, charge = FALSE,emass = 0.00054858, plotit = TRUE, algo=1, rel_to = 0, verbose = TRUE,return_iso_calc_amount = FALSE)

envelope(pattern=patt, ppm = FALSE, dmz = "get", frac = 1/4, env = "Gaussian",resolution = 5e+03, plotit = TRUE, verbose = TRUE)


library(Rdisop)
Mm <- getMolecule("C60H30O30", z=0, elements=initializeElements(c("C","H","O")))
getMass(Mm)
getFormula(Mm)
getIsotope(Mm)
getIsotope(Mm, index=1)
getValid(Mm)
isotopes <- getIsotope(Mm, seq(1,10))
isotopes
plot(isotopes[1,], isotopes[2,],type = "h",xlab = "m/z",ylab = "abundance",main = "C60H30O30",lwd=2)
abline(h=0,lty=2)


library(BRAIN) # bioc
aC = list(C=60, H=30, N=0, O=30, S=0)
nrPeaks = calculateNrPeaks(aC)
nrPeaks
res <- calculateIsotopicProbabilities(aC = aC, stopOption="nrPeaks", nrPeaks = nrPeaks)
res


library(IsoSpecR)
IsoSpecify(molecule = c(C=10,H=22,O=1), stopCondition = .9999 )

```


Izotopové poměry

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(CIAAWconsensus)

## Normalizace údajů o izotopech platiny na platinu-195
normalize.ratios(platinum.data, "platinum", "195Pt")

## Konsenzuální poměry množství izotopů platiny 
df=normalize.ratios(platinum.data, "platinum", "195Pt") 
mmm(df$R, df$u.R)

## Izotopové poměry zinku z izotopových zastoupení 
x = c(0.48630, 0.27900, 0.04100, 0.18750, 0.00620) 
ux = c(0.00091, 0.00076, 0.00031, 0.00135, 0.00010) 
z = abundances2ratios(x,ux,ref=2); z
at.weight(z$R,z$R.cov,"zinc","66Zn")

## Atomová hmotnost a izotopové zastoupení iridia, které odpovídají poměru izotopů 191Ir/193Ir = 0,59471(13) 
at.weight(0.59471, matrix(0.00013^2), "iridium", "193Ir")

## Atomová hmotnost a izotopové zastoupení křemíku, které odpovídají izotopovým poměrům 28Si/29Si = 1,074(69) a 30Si/29Si = 260(11) s korelací 0,80 mezi oběma izotopovými poměry. 
ratios = c(1.074,260) 
r.cov = matrix(c(0.069^2,0.80*0.069*11,0.80*0.069*11,11^2),ncol=2,byrow=TRUE) 
at.weight(ratios, r.cov, "silicon", "29Si")

```


Vlastnosti vzduchu a vody

```{r}

# Slozeni atmosfery
library(marelac)
atmComp() 

## Hustota vzduchu 
library(marelac)
air_density(t = 25, P = 1.013253)
# t - temperature in °C, P - pressure in bars

## Hustota vody
library(rLakeAnalyzer)
water.density(wtr = 25, sal = 0)
# wtr - temperature in °C
library(CHNOSZ)
rho.IAPWS95(T = 298.15, P = 1, state="", trace=0)
# T - temperature in K, P - pressure in bars


library(IAPWS95)
library(marelac)
library(CHNOSZ)

library(CHNOSZ)
thermo()$OBIGT

```




Disociační konstanty

```{r echo=FALSE, message=FALSE, warning=FALSE}

#https://chem.libretexts.org/Ancillary_Materials/Reference/Reference_Tables/Equilibrium_Constants/E1%3A_Acid_Dissociation_Constants_at_25C


library(seacarb)
Kd1 = K1p(S=0, T=25, P=0, pHscale="T", kSWS2scale="x", warn="n")
pKa1 = -log10(Kd1); pKa1
Kd2 = K2p(S=0, T=25, P=0, pHscale="T", kSWS2scale="x", warn="n")
pKa2 = -log10(Kd2); pKa2
Kd3 = K3p(S=0, T=25, P=0, pHscale="T", kSWS2scale="x", warn="n")
pKa3 = -log10(Kd3); pKa3

library(seacarb)
bjerrum(K1=Kd1[1], K2=Kd2[1], K3=Kd3[1], phmin=1, phmax=13, by=0.1, conc=1,type="l", col=c(2,3,4,6), ylab="Relative concentration (%)", add=FALSE)


library(SolveSAPHE)
AK_PHOS_1_MILL95(t_k=298,s=0, p_bar=0)
AK_PHOS_2_MILL95(t_k=298,s=35, p_bar=0)
AK_PHOS_3_MILL95(t_k=298,s=0, p_bar=0)


library(AquaEnv)
data(PhysChemConst)
PhysChemConst

```


```{r}

library(seqinr)
data(aaindex)
aaindex

a("Glu")
aaa("A")


library(CHNOSZ)
thermo = thermo()
thermo$protein
pf <- protein.formula(thermo$protein); pf
iA <- info("alanine")
info(iA)$formula
as.chemical.formula(makeup(iA))

```



Převody jednotek

```{r echo=FALSE, message=FALSE, warning=FALSE}
library(measurements)
conv_unit_options # pouzivane jednotky

# length
conv_unit(1, from="mm", to="m")
conv_unit(1, from="nm", to="angstrom")

# pressure
conv_unit(101, from="Pa", to="mmHg")
conv_unit(1, from="atm", to="Pa")

# temperature
conv_unit(100, from="C", to="K")

# energy
conv_unit(500, from="cal", to="J")
conv_unit(1, from="erg", to="J")

#volume
conv_unit(100, from="l", to="m3")

# mass
conv_unit(100, from="g", to="mg")
conv_unit(100, from="ug", to="g")

# speed
conv_unit(1, from="km_per_hr", to="m_per_sec")
conv_unit(1, from="light", to="km_per_hr")
conv_unit(1, from="light", to="m_per_sec")


### multiunits
conv_multiunit(x = 1, from="ug / l", to="g / m3")
conv_multiunit(x = 1, from="mg / l", to="kg / m3")
conv_multiunit(x = 1, from="cal / kg", to="J / g")


## Prevod teplot
library(weathermetrics)
convert_temperature(temperature=0, old_metric="c", new_metric="k", round = 2)
# metric = "c","f","k"
celsius.to.fahrenheit(T.celsius, round = 2)
celsius.to.kelvin(T.celsius, round = 2)
fahrenheit.to.kelvin(T.fahrenheit, round = 2)
fahrenheit.to.celsius(T.fahrenheit, round = 2)
kelvin.to.celsius(T.kelvin, round = 2)
kelvin.to.fahrenheit(T.kelvin, round = 2)


#### nasobky jednotek SI
library(sitools)
f2si(10000)
f2si(0.023, unit="l")
f2si(3.5e-5, unit="mol")
numbers <- c(1e5, 3.5e-12, 0.004) 
f2si(numbers, unit="g")

sapply(numbers,function(x){f2si(x, unit="g")})


# Prevod casovych jednotek

library(datetime)

xx = as.minute(60) # 60 min
as.second(xx)
as.hour(xx)
as.day(xx)

library(lubridate)
dur = duration(week=0,day=0.04166667,hours=1,minutes=0,seconds=3600)
as.numeric(dur, "minutes")


# Opticke jednotky
library(photobiology)
T2A(0.99, action = NULL, byref = FALSE, clean = TRUE) # transmitance na absorbanci
T2Afr(0.99, Rfr = 1/4, byref = FALSE, clean = TRUE) # transmitance na absorptanci (s danou mirou reflektance vzorku)
A2T(0.01, action = NULL, byref = FALSE, clean = TRUE) # absorbance na transmitanci


```


Plynný systém mění svůj objem o 1200 ml za konstantního vnějšího tlaku 30 atm. Jakou práci vykoná plyn při expanzi? Vyjádřete v různých energetických jednotkách.

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(measurements)

# W = p . dV
dV = conv_unit(1200, from="ml", to="m3"); dV
p = conv_unit(30, from="atm", to="Pa"); p
W = p*dV; W # [J]
W1 = conv_unit(W, from="J", to="cal"); W1 # [cal]
W2 = conv_unit(W, from="J", to="Wsec"); W2 # [W.s]
W3 = conv_unit(W, from="J", to="erg"); W3 # [erg]

```


Ocelová láhev o objemu 20.5 l je naplněna kyslíkem. Při 17 °C je tlak v láhvi 87 atm. Vypočítejte hmotnost kyslíku.

```{r}

library(marelac)
R = as.numeric(Constants$gasCt2[1]); R
Constants$gasCt2[2]
library(measurements)
P = conv_unit(87, from="atm", to="Pa"); P
Tk = conv_unit(17, from="C", to="K"); Tk
V = conv_unit(20.5, from="l", to="m3"); V
library(CHNOSZ)
M = mass("O2"); M # g/mol

m = P*V*M/(R*Tk)
m = ceiling(m); m # g
mk = conv_unit(m, from="g", to="kg")
mk = round(mk,2); mk # kg

paste("Hmotnost kyslíku je", m, "g, resp.", mk, "kg.")

```


Vypočítejte molekulovou hmotnost benzenu, váží-li 600 ml jeho par při 87 °C a tlaku 624 mm Hg 1.3 g.

```{r}

m = 1.3 # g
library(marelac)
R = as.numeric(Constants$gasCt2[1]); R
Constants$gasCt2[2]
library(measurements)
P = conv_unit(624, from="mmHg", to="Pa"); P
Tk = conv_unit(87, from="C", to="K"); Tk
V = conv_unit(600, from="ml", to="m3"); V

M = m*R*Tk/(P*V); M # g/mol

library(CHNOSZ)
Mb = mass("C6H6"); Mb # g/mol

paste("Zjištěná molární hmotnost benzenu je", round(M,0), "g/mol.")

```


Pokus trval celkem 7658 minut. Kolik je to dní, hodin a minut?

```{r echo=FALSE, message=FALSE, warning=FALSE}

tt = 7658 # [min]

library(lubridate)
tim = seconds_to_period(tt*60); tim
library(datetime)
tim = seconds_to_period(as.second(as.minute(tt))); tim

tn = regmatches(tim, gregexpr("[[:digit:]]+", tim))
tn = as.numeric(unlist(tn))
t1 = tn[1]
t2 = tn[2]
t3 = tn[3]
t4 = tn[4]
  
paste("Pokus trval", t1, "dní,", t2, "hodin,", t3, "minut a", t4, "sekund.")

```



Zaokrouhlování čísel na daný pocet desetinných míst

```{r echo=FALSE, message=FALSE, warning=FALSE}

round(1234.567,2)
round(123.456,digits=2)
ceiling(1234.567)
floor(1234.567)
trunc(1234.567)

library(guf)
round_something(123.456, decimals = 2)
round_something("123.456", decimals = 2)

library(PKNCA)
roundString(3141.59265, digits = 3, sci_range = 0, sci_sep = "e")
roundString(3141.59265, digits = 3, sci_range = 0, sci_sep = "x10^")
roundString(3141.59265, digits = 3, sci_range = Inf, sci_sep = "e")

# vektory
x <- c(sort(sample(1.765:20, 10)))
round(x,digits=2)
ceiling(x)
floor(x)
trunc(x)

library(ch)
Round(x, n=2)
Round2(x, n=2)


# matice

X = matrix(x,5,2,byrow = TRUE)
round(X,digits=2)
ceiling(X)
floor(X)
trunc(X)

Round(X, n=2)
Round2(X, n=2)

```


Nejistoty a chyby měření

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(errors)
xe <- set_errors(3.602, 0.008)
options(errors.notation = "plus-minus") 
print(xe, digits = 2)
options(errors.notation = "parenthesis")
print(xe, digits = 2)

# elementarni naboj
e <- set_errors(1.6021766208e-19, 0.0000000098e-19) 
options(errors.notation = "plus-minus") 
print(e, digits = 2)
options(errors.notation = "parenthesis")
print(e, digits = 3)

```


Ze zásilky kaprolaktamu bylo odebráno 10 vzorku a byl u nich stanoven bod tání. Vypočítejte průměrnou hodnotu bodu tání v zásilce a její směrodatnou odchylku.

```{r echo=FALSE, message=FALSE, warning=FALSE}

xc = c(68.5, 68.7, 68.3, 68.8, 68.5, 68.2, 68.6, 68.4, 68.2, 68.7)
mean(xc)
sd(xc) # sd
sd(xc)/sqrt(length(xc)) # str chyba prumeru

e <- set_errors(mean(xc), sd(xc)/sqrt(length(xc))) 
options(errors.notation = "plus-minus"); print(e, digits = 2)

```


Sireni chyb

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)

# z prumeru a nejistoty, bez uvedení poctu stupnu volnosti, resp. poctu opakovani

x <- c(5, 0.01) 
y <- c(1, 0.01) 

EXPR1 <- expression(x/y) 
DF1 <- cbind(x, y)
RES1 <- propagate(expr = EXPR1, data = DF1, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000) 
RES1
RES1$data
RES1$prop
RES1$sim
summary(RES1)
plot(RES1)

EXPR <- expression(a^b*x)
a = c(5, 0.1)
b = c(10, 0.1)
x = c(1, 0.1)
DAT <- cbind(a, b, x)
(res <- propagate(EXPR, DAT))
res$data
res$prop
res$sim
summary(res)
plot(res)


# z prumeru a nejistoty, s uvedením poctu stupnu volnosti
EXPR2 <- expression(x/y) 
x <- c(5, 0.01, 12) 
y <- c(1, 0.01, 5) 
DF2 <- cbind(x, y) 
RES2 <- propagate(expr = EXPR2, data = DF2, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES2
RES2$data
RES2$prop
RES2$sim
summary(RES2)
plot(RES2)


# Výpocet pomocí intervalu 

EXPR3 <- expression(C * sqrt((520 * H * P)/(M *(t + 460))))
H <- c(64, 65)
M <- c(16, 16.2)
P <- c(361, 365)
t <- c(165, 170)
C <- c(38.4, 38.5)
DAT3 <- makeDat(EXPR3)
interval(DAT3, EXPR3, seq = 2)

EXPR5 <- expression(x^2 - x + 1)
x <- c(-2, 1)
curve(x^2 - x + 1, -2, 1)
DAT5 <- makeDat(EXPR5)
interval(DAT5, EXPR5, seq = 2)


library(metRology)

data(GUM.H.1)
GUM.H.1

## a simple uncertainty analysis for the product of two quantities
GUM(c("x1","x2"),c(2.3,1.1),c(0.030,0.015),c(5,9999),"x1*x2")

## a simple uncertainty analysis for the product of two quantities
GUM.validate(c("x1","x2"), c(2.3,1.1), c(0.030,0.015), c(5,9999), c("A","B"),c("Normal","Rectangular"),"x1*x2")


expr <- expression(a+b*2+c*3+d/2)
x <- list(a=1, b=3, c=2, d=11)
u <- lapply(x, function(x) x/10) # nejistota 10 %
u.expr<-uncert(expr, x, u, method="NUM")
u.expr
# function method
f <- function(a,b,c,d) a+b*2+c*3+d/2
u.fun<-uncert(f, x, u, method="NUM")
u.fun
# formula method
u.form<-uncert(~a+b*2+c*3+d/2, x, u, method="NUM")
u.form


# s korelaci
u.cor<-diag(1,4)
u.cor[3,4]<-u.cor[4,3]<-0.5
u.cor
# num
u.formc<-uncert(~a+b*2+c*3+d/2, x, u, method="NUM", cor=u.cor)
u.formc
# Monte Carlo
u.formc.MC<-uncert(~a+b*2+c*3+d/2, x, u, method="MC", cor=u.cor, B=200)
u.formc.MC


expr <- expression(a/(b-c))
x <- list(a=1, b=3, c=2)
u <- lapply(x, function(x) x/20)
set.seed(403)
u.invexpr<-uncertMC(expr, x, u, distrib=rep("norm", 3), B=999, keep.x=TRUE )
u.invexpr

```


Vypočítejte hustotu a její chybu pro látku, u níž byla opakovaným měřením stanovena hmotnost 6.824 (0.008) g a objem 3.03 (0.01) ml.

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)

EXPR2 <- expression(m/V) 
m = c(6.824, 0.008) #[g] 
V = c(3.03, 0.01) #[ml]
DF2 <- cbind(m, V) 
RES2 <- propagate(expr = EXPR2, data = DF2, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES2
RES2$data
RES2$prop
RES2$sim
library(errors)
e <- set_errors(RES2$sim[1], RES2$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)

```


Na píst o průměru 200 (0.05) mm působí pára tlakem 8.2 (0.1) atm. Jakou silou působí pára na píst?

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)
library(measurements)

d = c(200, 0.05) # [mm] na [m]
d = as.vector(conv_unit(c(200, 0.05), from="mm", to="m"))
p = c(8.2, 0.1)  # [atm] na [Pa]
p = as.vector(conv_unit(c(8.2, 0.1), from="atm", to="Pa"))
pi
EXPR7 <- expression(p*3.141593*(d/2)^2)
DF7 <- cbind(d, p) 
RES7 <- propagate(expr = EXPR7, data = DF7, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES7
RES7$data
RES7$prop
RES7$sim
library(errors)
e <- set_errors(RES7$sim[1], RES7$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)

```


Vypočítejte koeficient viskozity roztoku glycerinu Stokesovou metodou

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)

r = c(0.0112, 0.0001) # polomer kulicky [cm]
l = c(31.23, 0.05) # draha kulicky za cas t [cm]
t = c(62.1, 0.2) # cas [s]
g = c(980.1,0) # tihove zrychleni [cm/s2] z tabulek
d0 = c(13.55,0) # hustota kulicky [g/cm3] z tabulek
d = c(1.28,0) # hustota roztoku [g/cm3] z tabulek
EXPR3 <- expression((2/9)*g*(((d0-d)*r^2)/l)*t)
DF3 <- cbind(r, l, t, g, d0, d) 
RES3 <- propagate(expr = EXPR3, data = DF3, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES3
RES3$data
RES3$prop
RES3$sim
EXPR4 <- expression((2/9)*980.1*(((13.55-1.28)*r^2)/l)*t)
DF4 <- cbind(r, l, t) 
RES4 <- propagate(expr = EXPR4, data = DF4, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES4$data
RES4$prop
RES4$sim
library(errors)
e <- set_errors(RES4$sim[1], RES4$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)

```


Vypočítejte koeficient viskozity roztoku glycerinu pomocí kapilárního viskozimetru.

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)
library(measurements)

p = c(20.12, 0.01) # tlak na vytoku z kapilary [mm Hg] to [Pa]
p = as.vector(conv_unit(c(20.12, 0.01), from="mmHg", to="Pa"))
r = c(0.570, 0.003) # polomer kapilary [mm] to [m]
r = conv_unit(c(0.570, 0.003), from="mm", to="m")
l = c(10.526, 0.005) # delka kapilary [mm] to [m]
l = conv_unit(c(10.526, 0.005), from="mm", to="m")
V = c(5.025, 0.001)# objem kapaliny vytekle za cas t [cm3] to [m3]
V = conv_unit(c(5.025, 0.001), from="cm3", to="m3")
t = c(27.34, 0.02) # cas [s]
pi
EXPR5 <- expression((3.141593*p*(r^4)*t)/(8*V*l))
DF5 <- cbind(p, r, l, V, t) 
RES5 <- propagate(expr = EXPR5, data = DF5, type = "stat", do.sim = TRUE, verbose = TRUE, nsim = 100000)
RES5 # [kg / m.s]
RES5$data
RES5$prop
RES5$sim
# conv_multiunit(x = RES5$sim[1:2], from="kg / m", to="g / cm")
library(errors)
e <- set_errors(RES5$sim[1], RES5$sim[2]) 
options(errors.notation = "plus-minus"); print(e, digits = 1)
options(errors.notation = "parenthesis"); print(e, digits = 1)

```


Součin rozpustnosti stříbrné soli AgX má hodnotu Ks = 4.0 (0.4) x 10^-8. Jaká je chyba vypočtené rovnovážné koncentrace stříbrných iontů ve vode?

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(propagate)
EXPR <- expression(x^0.5)
x = c(4.0, 0.4)
DAT <- data.frame(x)
res <- propagate(EXPR,DAT)
# NEFUNGUJE pro jednu promennou.


library(metRology)

GUM("Ks",4.0e-8,0.4e-8,1,"sqrt(Ks)",sig.digits.U = 2) 

expr <- expression(a^0.5)
x <- list(a=4.0e-8)
u <- lapply(x, function(x) x/10) # nejistota 10 %
u.expr<-uncert(expr, x, u, method="NUM")
u.expr
# function method
f <- function(a) a^0.5
u.fun<-uncert(f, x, u, method="NUM")
u.fun
# formula method
u.form<-uncert(~a^0.5, x, u, method="NUM")
u.form

```


Nejistoty a korelace

```{r echo=FALSE, message=FALSE, warning=FALSE}

library(errors)

x = c(0.9719378, 1.9840006, 2.9961830, 4.0123346, 5.0012799)
errors(x) = c(0.01, 0.01, 0.01, 0.01, 0.01)
y = c(0.9748992, 1.9627805, 2.9935831, 3.9921237, 4.9612555)
errors(y) = c(0.02, 0.02, 0.02, 0.02, 0.02)

# bez korelace
correl(x, y) =  c(0, 0, 0, 0, 0) 
z <- x / y; z

# s korelací
correl(x, y) =  c(0.8864282, 0.9761841, 0.9140209, 0.9496266, 0.9911837)
z_correlp <- x / y; z_correlp

correl(x, y) =  c(-0.8864282, -0.9761841, -0.9140209, -0.9496266, -0.9911837)
z_correln <- x / y; z_correln

```