# Evgenii Rudnyi
# requires matplotlib (it includes numpy)
# python -m pip install -U matplotlib
# http://blog.rudnyi.ru/ru/2024/02/reversible-heat-exchange.html

import numpy
import pylab
import math
import numbers

def PlotResultTogether(X, Y, 
		FileName = 'plot', 
		CommonTitle = '', 
		FunctionX = lambda x: x, 
		FunctionY = lambda x: x,
		LogX = False,
		LogY = False,
		LimX = False,
		LimY = False,
		Options = False,
		Grid = False,
		LabelX = False,
		LabelY = False,
		Legend = False,
		Location = 0):
	"""makes plots in the same figure. 
	
	Options: FileName for the base name, CommonTitle for a title,
	FunctionX and FunctionY to modify series values, LogX and LogY to
	turn on a logarithmic scale for X and Y axes, Options is a list of
	dictionaries with matplotlib properties for each series."""

	pylab.title(CommonTitle)
	ax = pylab.gca()
	if LogX:
		ax.set_xscale("log")
	if LogY:
		ax.set_yscale("log")
	pylab.grid(Grid)
	if LimX:
		pylab.xlim(LimX)
	if LimY:
		pylab.ylim(LimY)
	if LabelX:
		pylab.xlabel(LabelX)
	if LabelY:
		pylab.ylabel(LabelY)
	kk = 0
	for i in range(len(Y)):
		if isinstance(X[0], numbers.Number):
			x = [FunctionX(k) for k in X]
		else:
			x = [FunctionX(k) for k in X[i]]
		y = [FunctionY(k) for k in Y[i]]
		if Options == False:
			pylab.plot(x, y)
		else:
			pylab.plot(x, y, **Options[kk])
			kk = kk + 1
			if kk == len(Options):
				kk = 0
	if Legend != False:
		pylab.legend(Legend, loc = Location)
	pylab.savefig(FileName, dpi=125)
	pylab.clf()

def ComputeChangeS(x,T_ini):
	S = 0
	N = len(x)
	for i in range(N):
		S = S + math.log((x[i] + 273.15)/(T_ini+273.15))/N
	return S

def run(N):
	n = numpy.full((N,),100.)
	m = numpy.zeros((N,))
	for i in range(N):
		for j in range(i + 1):
			val = (n[-i+j-1] + m[j])/2
			m[j] = n[-i+j-1] = val
	for i in range(N - 1):
		for j in range(N - 1 - i):
			val = (n[j] + m[j + i + 1])/2
			m[j + i + 1] = n[j] = val
	delS1 = ComputeChangeS(n, 100)
	delS2 = ComputeChangeS(m, 0)
	return numpy.mean(n), numpy.mean(m), delS1 + delS2, n, m


arrN = [5, 10, 50, 100, 500, 1000, 5000, 10000]
delS = []
meann = []
meanm = []
vecn = []
vecm = []
for i in arrN:
	print(i)
	a, b, c, d, e = run(i)
	print(a, b, c)
	meann.append(a)
	meanm.append(b)
	delS.append(c)
	vecn.append(d)
	vecm.append(e)
	PlotResultTogether(range(i), [d, e], Legend = ['hot', 'cold'], FileName='temp_' + str(i))

ia = 0
ib = 4
ic = 7 
na = arrN[ia]
nb = arrN[ib]
nc = arrN[ic]
a = numpy.linspace(0., 1., na)
b = numpy.linspace(0., 1., nb)
c = numpy.linspace(0., 1., nc)
legend = ['hot_' + str(na), 'cold_'  + str(na), 'hot_' + str(nb), 'cold_'  + str(nb), 'hot_' + str(nc), 'cold_'  + str(nc)]
PlotResultTogether([a, a, b, b, c, c], [vecn[ia], vecm[ia], vecn[ib], vecm[ib], vecn[ic], vecm[ic]], Legend = legend, LabelX = 'Length', LabelY='Temperature', FileName='temp_together')

PlotResultTogether(arrN, [meann, meanm], Legend = ['hot', 'cold'], FileName='meanT', LogX = True)
PlotResultTogether(arrN, [delS], Legend = ['delS'], FileName='delS', LogX = True, LogY = True)