# -*- coding: utf-8 -*-
#TFY4165 Termisk fysikk 2022: Obligatorisk øving nr 1
#Simulering av loppebefengte hunder, se PCH kapittel 4.8.

import matplotlib.pyplot as plt
import numpy as np

#N = antall lopper
N = 400
#NT = antall tidssteg
NT = N*100
#c = eksponentiell decay-faktor
c = 1/N
#Nullstiller "tidstabeller" for antall lopper
#på hhv hund A og B
NA = np.zeros(NT)
NB = np.zeros(NT)
NAexp = np.zeros(NT)
tid = np.zeros(NT)
#Alle lopper på hund A i starten.
#Paa hund A: loppested = 1;  Paa hund B: loppested = 0
loppested = np.ones(N)
#print(loppested)
#Tabell med tilfeldig valgt loppe som hopper.
#NT tilfeldige heltall mellom 0 og N-1.
#Dvs, loppene er nummerert fra 0 til N-1,
#i god pythonånd.
hopper = np.random.randint(0,high=N,size=NT)
#print(N,NT)
#print(hopper)
#Alle på hund A i starten
NA[0]=N
NB[0]=0
#Analytisk NA(t)
NAexp[0] = (N/2)*(1+np.exp(-2*c*0))
tid[0] = 0
for i in range(1,NT,1):
    if loppested[hopper[i]] == 1:
        loppested[hopper[i]] = 0
        NB[i]=NB[i-1]+1
        NA[i]=NA[i-1]-1
    else:
        loppested[hopper[i]] = 1
        NB[i]=NB[i-1]-1
        NA[i]=NA[i-1]+1
    
    NAexp[i] = (N/2)*(1+np.exp(-2*c*i))
    tid[i] = i

#Plotting
plt.figure('Lopper')
#plt.title('Lopper')
plt.plot(tid,NA,tid,NAexp)
plt.ylabel('NA(t)/N')
plt.xlabel('t')
plt.xlim(0,NT)
plt.ylim(0,N)
plt.grid()
plt.show()