%%Conduction tutorial Summer 2014
%%Dimensionless flux qy vs. x at y=0

%Read data from Excel file
clear all; clc;

%ANSYS solution
x_ANSYS = xlsread('boundary_flux_tutorial.xlsx', 1, 'A4:A105');
qy_ANSYS = xlsread('boundary_flux_tutorial.xlsx', 1, 'B4:B105'); 
        %10x20 mesh
qy_ANSYS2 = xlsread('boundary_flux_tutorial.xlsx', 1, 'C4:C105'); 
        %20x40 mesh   bias factor = 20
qy_ANSYS3 = xlsread('boundary_flux_tutorial.xlsx', 1, 'D4:D105'); %10x20 mesh
        %40x80 mesh   bias factor = 20
%Theory: 100 term series solution
% x_theory = xlsread('boundary_flux.xlsx', 1, 'F3:F103');
% qy_theory = xlsread('boundary_flux.xlsx', 1, 'G3:G103');

figure(1); clf;
h=plot(x_ANSYS, qy_ANSYS,'-k');
set(gca,'Box','on','LineWidth',2,'FontName','Helvetica',...
	'FontSize',14);
set(h,'LineWidth',2);
hold on;
h=plot(x_ANSYS, qy_ANSYS2,'--r');
set(h,'LineWidth',2);
h=plot(x_ANSYS, qy_ANSYS3,'-.b');
set(h,'LineWidth',2);
xlabel('x^*'); ylabel('q^+');
axis square;
legend('ANSYS 10x20 mesh','ANSYS 20x40 mesh with bias',...
     'ANSYS 40x80 mesh with bias','Location','NorthWest');
title('Dimensionless Flux at Bottom Boundary');
axis([0 1.05 0 20]);

%Overall dimensionless flux crossing bottom surface
overall_flux = trapz(x_ANSYS, qy_ANSYS3)