plot_case.py

import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

case_name = 'caso1'

center_data = np.loadtxt(case_name, skiprows=4, max_rows=401, delimiter='\t')
surface_data = np.loadtxt(case_name, skiprows=408, max_rows=401, delimiter='\t')

Tmax = 523
Tmin = 278

Bi_f = 0.1
C = 10

center_data[:,1] = (center_data[:,1] - Tmin)/(Tmax - Tmin)
center_data[:,0] /= np.max(center_data[:,0])

surface_data[:,1] = (surface_data[:,1] - Tmin)/(Tmax - Tmin)
surface_data[:,0] /= np.max(surface_data[:,0])

ideal_case = np.zeros_like(surface_data)
ideal_case[:,0] = surface_data[:,0]
ideal_case[:,1] = (np.cosh(C*(ideal_case[:,0])*np.sqrt(Bi_f)) + np.sqrt(Bi_f)*np.sinh(C*(ideal_case[:,0])*np.sqrt(Bi_f))) / (np.cosh(C*np.sqrt(Bi_f)) + np.sqrt(Bi_f)*np.sinh(C*np.sqrt(Bi_f)))
    
fig = plt.figure()
ax = fig.add_subplot(111)
    
ax.plot(center_data[:,0], center_data[:,1], '-.r', label='Centre')
ax.plot(surface_data[:,0], surface_data[:,1], '--b', label='Surface')
ax.plot(ideal_case[:,0], ideal_case[:,1], '-k', label='1D solution')

ax.set_xlabel('$x*$')
ax.set_ylabel('$\\theta*$ ', rotation = 0)

ax.autoscale(tight=True)
ax.set_xticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])
ax.set_yticks([0.0, 0.2, 0.4, 0.6, 0.8, 1.0])

ax.grid()
ax.legend()
ax.set_title('Case 1: Straight fin (Bi = 0.1, C = 10)')

fig.savefig(case_name + '.png')