%/* program to evaluate exact sensitivity coefficients for
%   the constant property case.

%   the sensitivity coefficients are for the HTC, where the
%   boundary condition at x=0 is
%		q = HTC*(T_inf - T(0))
%
%
%*/

% ported from KAW original .cpp code from "CastingProject\C++\exact_sens"
%   conversion June 2010
function [time, fvals, sensvals] = compute_sens(h,L,k,C,dt,tfinal,N_terms,nsensors,theta_i, x)
	Bi = h*L/k;
    global Bi_global;
	Bi_global = Bi;
	Cn = zeros(N_terms);
	zeta = get_eigenvalues(N_terms);
	Cn = compute_Cn(N_terms,zeta);
    [time, fvals, sensvals ] = evaluate(x,dt,tfinal,k,C,L,N_terms,zeta, Cn, theta_i, h);
end	


function [Cn] = compute_Cn( N_terms, zeta)
	for i = 1:N_terms
		Cn(i) = 4*sin(zeta(i))/(2.*zeta(i) + sin(2.*zeta(i)));
    end
end

function [time, func_val, sens_val] = evaluate(x,dt,tfinal,k,C,L,N_terms,zeta, Cn, theta_i, h)
    fprintf(1,' Results for h = %f ', h);
	for i = 1:length(x)
		fprintf(1,'       x = %f',x(i));
    end
	fprintf(1,'\n');
    j = 0;
	for t = 0.:dt:1.0001*tfinal
        j = j + 1;
		Fo = k/C*t/L/L;
		fprintf(1,'  %f      %f   ',t,Fo);
		for i = 1:length(x)
			func = thetastar( 1. - x(i)/L , Fo,N_terms, zeta, Cn );
			sensitivity = d_thetastar_d_Bi( 1. - x(i)/L , Fo, N_terms, zeta, Cn  );
			fprintf(1,'    %f    %f   ',func*theta_i,sensitivity*L/k*theta_i*h);
            time(j) = t;
            func_val(j,i) = func*theta_i;
            sens_val(j,i) = sensitivity*h*L/k*theta_i;
        end
    	fprintf(1,'\n');
    end
end

function value = thetastar( xstar, Fo, N_terms, zeta, Cn )
    value = 0;
	for i = 1:N_terms
		value = value + Cn(i)*exp(-zeta(i)*zeta(i)*Fo)*cos(zeta(i)*xstar);
    end
end

function value = d_zeta_d_Bi( zeta_n )
	value = 0;
	value = 1./( tan(zeta_n) + zeta_n*sec( zeta_n )*sec( zeta_n ) );
end

function value = d_Cn_d_Bi( zeta_n )
	num = 4.*cos(zeta_n)*(2*zeta_n + sin(2.*zeta_n)) ...
			- 8.*sin(zeta_n)*(1. + cos(2.*zeta_n));
	denom = (2.*zeta_n + sin(2.*zeta_n))*(2.*zeta_n + sin(2.*zeta_n));
	value = num*d_zeta_d_Bi( zeta_n )/denom;
end

function value = d_thetastar_d_Bi( xstar, Fo, N_terms, zeta, Cn )
	value=0.;
	for i = 1:N_terms
		part(1) = d_Cn_d_Bi( zeta(i) )*cos(zeta(i)*xstar);
		part(2) = - 2.*Cn(i)*zeta(i)*Fo*d_zeta_d_Bi( zeta(i) )*cos(zeta(i)*xstar);
		part(3) = - Cn(i) * xstar *sin( zeta(i)*xstar )*d_zeta_d_Bi( zeta(i) );
		value = value + (part(1)+part(2)+part(3))* exp( - zeta(i)*zeta(i)*Fo );
    end
end

function value =  eigenfunction( x )
    global Bi_global;
	value = 0;
	value = x*tan(x) - Bi_global;
end

function zeta = get_eigenvalues(N_terms)
    options = optimset('TolX', 1.e-6 );
	for i = 1:N_terms
		xlo = (i-1)*pi;
	    xhi = (2*i-1)*pi/2.*.99999999;
%		fprintf(1,'  %d    %f    %f    %f    %f \n', i, xlo, xhi, eigenfunction(xlo), eigenfunction(xhi) );
		zeta(i) = fzero( @eigenfunction, [xlo xhi], options);
		fprintf(1,'  %d    %f  \n', i, zeta(i)/pi);
    end
end
