require("ggplot2")

require("xtable")

require("gridExtra")



require("scales")

require("pracma")



require("readxl")







stdev <- function(x) {

  s = sqrt( sum( (x-mean(x) )^2 )/length(x) )

  return(s)

}





#indirizzo

#C:/Users/Guglielmo/Desktop/università/ottica laboratorio/esperienza 1



#STABILITÃ  DEL LASER



dati_rumore <- read_excel("C:/Users/Guglielmo/Desktop/università/ottica laboratorio/esperienza 1/01b_Dati_Rumore_di_fondo.xlsx")

dati_rumore <- as.matrix.data.frame(dati_rumore)



dati_fluttuazione <- read_excel("C:/Users/Guglielmo/Desktop/università/ottica laboratorio/esperienza 1/01_Dati_Fluttuazioni_laser.xlsx")

dati_fluttuazione <- as.matrix.data.frame(dati_fluttuazione)



tempo_fluttuazione <- c(dati_fluttuazione[,1])

fluttuazione <- c(dati_fluttuazione[,2])

tempo_fluttuazione=tempo_fluttuazione[-c(1,2)] #tolgo i primi due elementi del vettore ridefinendolo

tempo_fluttuazione=as.numeric(tempo_fluttuazione) #conversione char --> numero

fluttuazione = fluttuazione[-c(1,2)]

fluttuazione = as.numeric(fluttuazione) #stessa cosa



#fare grafico

plot(tempo_fluttuazione, fluttuazione,

     cex = 0.5,

     main = "Fluttuazioni nel tempo")



#istogramma

hist(fluttuazione, xlim=a, main="Istogramma delle fluttuazioni", xlab="fluttuazione [mV]")



#dati sul rumore giorno/notte

tempo_rumore <- c(dati_rumore [,1])

rumore_accese  <- c(dati_rumore [,2])

rumore_spente  <- c(dati_rumore [,3])

tempo_rumore=tempo_rumore[-c(1,2)] #tolgo i primi due elementi del vettore

tempo_rumore=as.numeric(tempo_rumore) #cambio tipo char --> numero

rumore_accese=rumore_accese[-c(1,2)]

rumore_accese=as.numeric(rumore_accese)

rumore_spente=rumore_spente[-c(1,2)]

rumore_spente=as.numeric(rumore_spente)





#grafici

plot(tempo_rumore, rumore_accese,

     cex = 0.5, col = "red",

     main = "Rumore nel tempo")

plot(tempo_rumore, rumore_spente,

     cex = 0.5, col = "blue",

     main = "Rumore nel tempo")



#istogramma

hist(rumore_accese, main="Istogramma rumore luci accese", xlab="fluttuazioni [mV]")

hist(rumore_spente, main="Istogramma rumore luci spente", xlab="fluttuazione [mV]")





#LEGGE DI MALUS





pol2 <- read_excel("C:/Users/Guglielmo/Desktop/università/ottica laboratorio/esperienza 1/02_2Polaroid_cifre.xlsx")

pol2 <- as.matrix.data.frame(pol2)



#qui i dati sono già nuemeri double senza bisogno di conversioni

#inoltre i vettori sono già della corretta lunghezza

I <- 1/6 *(pol2[,2] + pol2[,3] + pol2[,4] + pol2[,5] + pol2[,6] + pol2[,7])



incertezzaI <- c(1:length(pol2[,1]))



for(i in 1:length(pol2[,1])){

  I

  incertezzaI[i] <- sqrt(   ( (pol2[i,2] - I[i])^2 + (pol2[i,3] - I[i])^2 +

                              (pol2[i,4] - I[i])^2 + (pol2[i,5] - I[i])^2 +

                              (pol2[i,6] - I[i])^2 + (pol2[i,7] - I[i])^2  ) /6   )

}



#sommo in quadratura incertezza con quella strumentale (ultimo digit)

for(i in 1:length(incertezzaI)){

  incertezzaI[i]=sqrt((incertezzaI[i])^2+(0.001)^2 + (0.02*I[i])^2)

}



#aggiungo l'offset

for(i in 1:length(incertezzaI)){

  I[i]=I[i]-0.032

  incertezzaI[i]=sqrt((incertezzaI[i])^2+(0.002)^2)#ultimo pezzo da togliere



}









Malus <- function(x, I0, phi, c){

  x <- x*pi/180 #così ho l'angolo in radianti

  I = I0 * cos(x+phi)^2 + c

  return(I)

}





pesi=1/(incertezzaI^2)

d <- data.frame(x=pol2[,1], y=I)

fit <- nls(y ~ Malus(x, I0, phi, c), data = d,

           start = list(I0 = 5, phi=pi/2, c=0), weights = pesi)  #valori più o meno vicini a quelli giusti

summary(fit)



I0=

phi=

c=



plot(pol2[,1], I,

     xlab = expression(theta), ylab = "Intensità [V]",

     main = "Legge di Malus con 2 Polaroid",

     pch = 16,

     cex = 0.6

     )

segments(pol2[,1], I+incertezzaI, pol2[,1], I-incertezzaI)

arrows(pol2[,1], I+incertezzaI, pol2[,1], I-incertezzaI, length=0.02, angle = 90, code = 3 )

curve(Malus(x, I0, phi,c), add=TRUE, col ='blue')





pol3 <- read_excel("C:/Users/Guglielmo/Desktop/università/ottica laboratorio/esperienza 1/03_3Polaroid_cifre.xlsx")

pol3 <- as.matrix.data.frame(pol3)





x=pol2[,1]

#chi2

chi2=0

for(i in 1:length(I)){

  chi2=chi2+(((I[i]-Malus(x,I0,phi,c)[i])^2)/(incertezzaI[i]^2))



}











#ricalcolo incertezze sulle x nel fit (metodo iterativo)-----------



x=pol2[,1]

sigmax=c(1:length(x))

for(i in 1:length(x)){

  sigmax[i]=1   #ci vorrebbe /sqrt(3)

}





cor_sigma=I0*(pi/180)*2*cos((pi/180)*x+phi)*sin((pi/180)*x+phi)*sigmax

incertezzaI_nuova=sqrt((incertezzaI)^2+(cor_sigma)^2)



#rifaccio fit





pesi=1/(incertezzaI_nuova^2)

d <- data.frame(x=pol2[,1], y=I)

fit <- nls(y ~ Malus(x, I0, phi,c), data = d,

           start = list(I0 = 5, phi=pi/2,c=0), weights = pesi)  #valori più o meno vicini a quelli giusti

summary(fit)



I0=

phi=

c=



plot(pol2[,1], I,

     xlab = expression(theta), ylab = "Intensità [V]",

     main = "Legge di Malus con 2 Polaroid",

     pch = 16,

     cex = 0.6

)

segments(pol2[,1], I+incertezzaI_nuova, pol2[,1], I-incertezzaI_nuova)

arrows(pol2[,1], I+incertezzaI_nuova, pol2[,1], I-incertezzaI_nuova, length=0.02, angle = 90, code = 3 )

curve(Malus(x, I0, phi,c), add=TRUE, col ='blue')