mec2 <- 510.999
r <- 0.05
ee <- 1.60e-19
c <- 299792348
eV <- 6.24e18

table <- read.table('dades/dades.txt', sep="\t", dec='.', header=TRUE)
table$LogIntensitatReal <- log(table$IntensitatReal)
table$CampMagneticSq <- table$CampMagnetic^2
table$E <- (sqrt((ee*c*(table$CampMagnetic/1000)*r*eV/1000)^2 + mec2^2) - mec2)/1000
table$ESq <- table$E^2

tablefit <- table[8:17,]

fit <- lm(IntensitatReal ~ E + ESq, data = tablefit)
summary(fit)
fit_coef <- coef(fit)

svg(file='output/p1b.svg', width=6, height=5)
plot(table$E, table$IntensitatReal, xlab="Energia (MeV)", ylab="Intensitat (s^-1)")
curve(fit_coef[1] + fit_coef[2]*x + fit_coef[3]*x*x, col = "mediumorchid2", add = TRUE)
dev.off()

svg(file='output/p1b2.svg', width=6, height=5)
plot(table$CampMagnetic, table$E, xlab="Camp magnètic (mT)", ylab="Energia (MeV)")
dev.off()

Emax <- -fit_coef[2]/(2*fit_coef[3])
Emax