# PYTHON & PASCO: Newtons avkjølingsmetode
# Versjon: 15.6.2021
# Kontakt: koding@astronomen.no
from pylab import *

# -------------------- Verdier for simulering og regresjon ---------------------
# k har mye å si for simuleringen
# temp_rom hvor fort temperaturendringen skjer.
k = 0.01              # Avkjølingskonstant i s^(-1)
temp_rom = 23.0       # Temperatur i rommet i grader C.

# Modell for temperaturendring.
def temp_endring(temp):          # Funksjon som regner temperaturendringen
    endring = -k*(temp-temp_rom) # Endring proporsjonal med temp-forskjell
    return endring               # Returnerer endringen.

# --------------------------- Importering av data ------------------------------
# Lister for lagring av data
tid_data = []
temp_data = []
temp_reg_data = []
# Innlesing av data
with open('avkjoling.csv') as file:       # Åpner datafilen fra PASCO
    next(file)      # Hopper over overskrift-linjen
    # Løkke over hver linje i datafilen. Linjer leses som tekst/streng.
    for line in file:
        # Klargjør data for lagring i lister
        line = line.replace(",",".")      # Bytter desimaltegn
        line = line.replace("\n","")      # Fjerner uønskede linjeskift
        line = line.split(";")            # Splitter streng til liste
        # Tallene i lista står som tekst, konverterer til tall og lagrer
        tid_data.append(float(line[1]))   # Tidsdata lagres i egen liste
        temp_data.append(float(line[2]))  # Temperaturdata lagres i liste
        temp_reg_data.append(float(line[2]) - temp_rom)

# --------------------------- Regresjon på data --------------------------------
# Kurvetilpasning med polyfit(x,y,polynomgrad)
# Vi bruker at en eksponensiell modell er lineær med logaritmisk skala.
# Vi korrigerer for romtemperaturen slik at regresjonen fungerer
logaritmisk_temp = log(temp_reg_data)       # Andreaksen gjøres logaritmisk
a,b = polyfit(tid_data,logaritmisk_temp,1)  # Lineær regresjon på log-data
temp_reg = []                             # Liste for y-verdiene fra regresjon
for t_verdi in tid_data:                  # Løkke over tidspunkter i data
    temp = exp(a*t_verdi + b)             # Kalkulerer temperatur i et tidspunkt
    temp_reg.append(temp + temp_rom)      # Lagrer temperatur med romtemperatur

# -------------------------- Simulering av data --------------------------------
# Startverdier
tid = tid_data[0]       # Starttid settes til første tidsverdi i data
temp = temp_data[0]     # Starttemperatur settes til første temperatur i data

# Simuleringsteknisk
dt = 0.001            # Lengde på tidssteg, s
tid_slutt = tid_data[-1]    # Sluttid for simulering er siste tidsdata
tid_sim = [tid]       # Liste for lagring av tid
temp_sim = [temp]     # Liste for lagring av temperatur

# Løkke som beregner ny temperatur og tid
while tid < tid_slutt:
    temp = temp + temp_endring(temp)*dt
    tid = tid + dt

    temp_sim.append(temp)
    tid_sim.append(tid)
# ------------------------ Plotting av resultater ------------------------------
# Data fra PASCO plottes som punkter, regresjon og modell som grafer
plot(tid_data,temp_data,".")    # Plotter data som punkter
plot(tid_data,temp_reg)         # Plotter regresjonsgraf
plot(tid_sim,temp_sim)          # Plotter simulerte data
title("Newtons avkjølingslov, data, regresjon og simulering")
xlabel("Tid / s")               # x-akse tittel
ylabel("Temperatur / Celcius")  # y-akse tittel
legend(["Data","Regresjon","Simulering"]) # navngivning av punkter og grafer
grid()   # Tegner rutenett
show()   # Vis plottet