#/usr/bin/python2.4
"""Il Pianeta delle Margherite."""


### Parametri del pianeta ###

Es=100.		# Energia in ingresso per unita' di superficie
Var=90; Vel=1	# Variazione dell'energia e velocita' di variazione
C=10.**11	# Capacita' termica del pianeta
CSB=10**-8	# Costante di Stefan-Boltzman
S=10.**10	# Superficie totale
As=0.5		# Assorbanza per la superficie vuota


### Parametri delle margherite ###

An=0.999	# Assorbanza per le margherite nere
Ab=0.001	# Assorbanza per le margherite bianche
Tb=303		# Temperatura ideale per le margherite bianche
Tn=293		# Tempertatura ideale per le margherite nere
Rb=20		# Tolleranza delle margherite bianche
Rn=20		# Tolleranza delle margherite nere
d=10.**-2	# Mortalita' per iterazione
f=1.4		# Fattore di crescita 
e=1.1		# Fattore di fertilita' esponenziale


### Condizioni iniziali ###

T=250.		# Temperatura iniziale del pianeta
B=500.		# Superficie iniziale di bianchi
N=500.		# Superficie iniziale di neri


### Controlli programma ###

iter=3		# Numero di sequenze
niter=1000	# Iterazioni per sequenza


###
### Inizio programma ######################################
###


from math import sin, fabs,pi
from commands import getstatusoutput
from pylab import *

###
### Definisco le funzioni #################################
###

def proliferazione(x0,Ti,Ri):
	""" Controlla il tasso di riproduzione o di mortalita' delle margherite in funzione della temperatura del pianeta """
	
	confort=e**( -((fabs(T-Ti)))/(Ri))
	unconfort=e**( -2.3 )
	if confort>unconfort:
		x1=x0*( f*confort )
	else:
		x1=x0*(confort**0.3)
		
	x1=x1*(1-d)
	#if x1<1:
	#	x1=1
	return x1

def deltat():
	"""Calcola la variazione di temperatura in funzione della composizione superficiale del pianeta, tenendo conto dell'energia incidente e di quella irradiata """
	
	En = An*E*N
	Eb = Ab*E*B
	Esup = As*E*(S-B-N)
	Ein = En+Eb+Esup
	Is = CSB*(1-As)*(S-B-N)
	In = CSB*(1-An)*N
	Ib = CSB*(1-Ab)*B
	Eout = (Is+In+Ib)*(T**4)
	Dt = (Ein-Eout)/C
	return Dt

def deltat2():
	""" Calcola la variazione di temperatura per un pianeta di controllo senza margherite """
	
	Ein = As*E
	Eout = CSB*(1-As)*(T2**4)
	Dt=(Ein-Eout)*(S/C)
	return Dt

def export(fn):
	""" Riporta in grafico l'andamento delle popolazioni, delle temperature e dell'energia """
	
	source='marghe'+str(fn)+'.log'
	output=source.replace('log','png')
	data=load(source)
	T1=data[:,0]
	N=data[:,1]
	B=data[:,2]
	E=data[:,3]
	T2=data[:,4]
	figure(fn)
	subplot(212)
	plot(N,color='k')
	plot(B,color=0.6)
	xlabel('Iterazioni')
	ylabel('N e B')
	subplot(211)
	title('File '+source)
	plot(T1,color='b')
	plot(E,color='r')
	xlabel('Iterazioni')
	ylabel('T1, T2 ed E')
	plot(T2,color='y')
	savefig(output)
	print 'Grafico dall\'iterazione '+str((fn)*niter+1)+' alla '+str((fn+1)*niter)+ ' disegnato in '+output+'...'

###
### Inizializzo le variabili ####################################
###

T2=T
E=Es
i=0
fn=0
count=0
unconfort=e**( -2.3 )
if f*unconfort<1.001:
	f=1.001/unconfort

# Faccio pulizia dei file di log e dei grafici preesistenti
getstatusoutput('rm marghe*.log')
getstatusoutput('rm marghe*.png')

# Apro il primo log
log=open('marghe0.log','w')

###
### Loop principale ##############################################
###
while [True]:
	
	# Calcolo le nuove temperature
	T = T + deltat()
	T2 = T2 + deltat2()
	
	# Calcolo le nuove popolazioni di margherite
	Np = proliferazione(N,Tn,Rn)
	Bp = proliferazione(B,Tb,Rb)
	
	# Nel caso la superficie occupata ecceda quella esistente, ridistribuisco le popolazioni proporzionalmente
	if Bp+Np>S:
		N=S*Np/(Bp+Np)
		B=S*Bp/(Bp+Np)
	else:
		N=Np
		B=Bp
	
	
	i+=1
	count+=1
	
	# Scrivo i risultati sul file di log
	if count<niter:
		log.write(str(T)+'\t'+str(N)+'\t'+str(B)+'\t'+str(E)+'\t'+str(T2)+'\n')
	else:
		# Se la sequenza Ú terminata: chiudo il log, traccio il grafico, e apro un nuovo log
		log.close()
		export(fn)
		fn+=1
		if fn==iter:
			break
		log=open('marghe'+str(fn)+'.log','w')
		count=0
	
	# Faccio variare la potenza della stella
	E=Es+Var*sin(Vel*i*pi/180)

print "Simulazione conclusa."

getstatusoutput('eog marghe0.png')