% fdX12B10T0 function
% Revision History
% 1 2 2013 James V. Beck, Matthew P. Lempke and Filippo de Monte
% INPUTS:
% xd: dimensionless location starting at xd=0 and ending at xd=1
% td: dimensionless time starting at td=0
% A: desired accuracy (1E-A = 10^-A); A = 2, 3, ..., 15
% OUPUTS:
% Td: dimensionless temperature calculated at (xd,td) to desired accuracy A
% qd: dimensionless heat flux calculated at (xd,td) to desired accuracy A
% Calling Sequence:
% None called
function [Td,qd]=fdX12B10T0(xd,td,A)
if( nargin < 3 )
    A = 15;  % default accuracy
end
lengthx=length(xd);
lengtht=length(td);
Td=zeros(lengtht,lengthx); % Preallocating Arrays for speed
qd=zeros(lengtht,lengthx); % Preallocating Arrays for speed
for it=1:lengtht % Begin time loop
    td_it=td(it); % Set current time
    for ix=1:lengthx % Begin space loop
        xd_ix=xd(ix); % Set current space
        td_dev1=(1/(10*A))*(2-xd_ix)^2; % First deviation time
        td_dev2=(1/(10*A))*(2+xd_ix)^2; % Second deviation time
        if td_it == 0 % For time t=0 condition
            Td(it,ix)=0; % Set inital temperature
            qd(it,ix)=0; % Set inital heat flux
        elseif td_it <= td_dev1
        % Solution for first small times:
            Td(it,ix)=erfc(xd_ix/sqrt(4*td_it));
            qd(it,ix)=1/(sqrt(pi*td_it))*exp(-xd_ix^2/(4*td_it));
        elseif td_it > td_dev1 && td_it <= td_dev2
        % Solution for second small times:
            Td(it,ix)=erfc(xd_ix/sqrt(4*td_it))+...
            erfc((2-xd_ix)/sqrt(4*td_it));
            qd(it,ix)=1/(sqrt(pi*td_it))*exp(-xd_ix^2/(4*td_it))-...
            1/(sqrt(pi*td_it))*exp(-(2-xd_ix)^2/(4*td_it));
        else
            m_max=ceil(1/2+sqrt(A*log(10)/(td_it*pi^2))); %Set max.No.terms
            % Start X12B10T0 case for large times:
            Td(it,ix)=1; % steady-state temperature solution
            qd(it,ix)=0; % steady-state heat flux solution
            for m=1:m_max % Continue X12B10T0 case for large times
            % Series solutions:
                betam=(m-1/2)*pi; % Define eigenvalues
                Td(it,ix)=Td(it,ix)-2*exp(-betam^2*td_it)*...
                sin(betam*xd_ix)/betam;
                qd(it,ix)=qd(it,ix)+2*exp(-betam^2*td_it)*cos(betam*xd_ix);
            end % for m
        end % if td_it
        if xd_ix == 1 % For location x=1 condition
            qd(it,ix)=0; % Set x=1 boundary heat flux
        end % if xd_ix
    end % for ix
end 