%  Temperature profiles in 1D conduction 
function [Tspan, rspan] = Tprofile1D(varargin)
  plate=varargin{1};
  if length(varargin{2})==1
    qgen=varargin{2};
  else
    qgen=0;
  end
  r=plate{end};
  r1=r(1); r2=r(end);
  if length(r)>2 rspan=r;
  else rspan=linspace(r1,r2,101);
  if length(plate)<5
    T=varargin{end};
    [id, k]=plate{1:2};
    if id=='Rec' | id=='Bar'
      F=@(r,X)(X(1)*(r2-r)+X(2)*(r-r1))/(r2-r1);
      Tspan=F(rspan,T)+qgen/2/k*(F(rspan,[r1,r2].^2)-rspan.^2);
    elseif id=='Cyl'
      if r1==0 G=@(r,X)X(2);
      else G=@(r,X)(X(1)*log(r2./r)+X(2)*log(r./r1))/log(r2/r1);
      end
      Tspan=G(rspan,T)+qgen/4/k*(G(rspan,[r1,r2].^2)-rspan.^2);
    elseif id=='Sph'
      if r1==0 H=@(r,X)X(2)./r2;
      else H=@(r,X)(X(1)*(r2./r-1)+X(2)*(1-r1./r))/(r2-r1);
      end
      Tspan=H(rspan,[T(1)*r1,T(2)*r2])+qgen/6/k*(H(rspan,[r1,r2].^3)-rspan.^2);
    else
      error('wall id must = ''Rec'', ''Bar'', ''Cyl'' or ''Sph''');
    end
  else
    T=varargin{end-1};
    Tinf=varargin{end};
    [id, k, Acond, Pconv, h]=plate{1:5};
    if length(h)~=length(Tinf)
      error('Dimensions of h and Tinf are not equal!');
    end
    havg=sum(h)/length(h);
    Teff=sum(h.*Tinf)/sum(h);
    if id=='Num'
      l=Acond(r(end))/Pconv(r(end));
    else
      l=Acond/Pconv;
    end
    Bi=havg*l/k;
    if Bi > 0.1  warning = 'Tprofile1D: Bi > 0.1'
    end
    Tp=qgen*l^2/k/Bi;
    Theta=T(1:2)-Teff-Tp;
    if id=='Rec' | id=='Bar'
      L=abs(r2-r1);
      lam=sqrt(Bi).*L./l;
      lespan=lam*(rspan-r1)/L;
      Th=Theta(1)*cosh(lespan)+ ...
        (Theta(2)-Theta(1)*cosh(lam))/sinh(lam)*sinh(lespan);
      Tspan=Th+Tp+Teff;
     elseif id=='Cyl'
      lam=sqrt(Bi).*[r1,r2]./l;
      lam(find(lam==0))=max(lam)*1e-6;
      E=besseli(0,lam(2)).*besselk(0,lam(1))-besseli(0,lam(1)).*besselk(0,lam(2));
      D1=Theta(1)*besseli(0,lam(2))-Theta(2)*besseli(0,lam(1));
      D2=Theta(2)*besselk(0,lam(1))-Theta(1)*besselk(0,lam(2));
      lrspan=sqrt(Bi)./l.*rspan;
      lrspan(find(lrspan==0))=max(lrspan)*1e-6;
      Th=(D1*besselk(0,lrspan)+D2*besseli(0,lrspan))/E;
      Tspan=Th+Tp+Teff;
    elseif id=='Sph'
      R=plate{6};
      lam=(1-cos([r1,r2]/R))/2;
      lam(find(lam==0))=max(lam)*1e-6;
      B=Bi*(R/l)^2;
      [~, Fx1]=Fsph(B,lam(1));
      [~, F1_x1]=Fsph(B,1-lam(1));
      [~, Fx2]=Fsph(B,lam(2));
      [~, F1_x2]=Fsph(B,1-lam(2));
      E=Fx1*F1_x2-Fx2*F1_x1;

      E1=Theta(1)*F1_x2-Theta(2)*F1_x1;
      E2=Theta(2)*Fx1-Theta(1)*Fx2;
      lspan=(1-cos(rspan/R))/2;
      lspan(find(lspan==0))=max(lspan)*1e-6;
      [~,Flam]=Fsph(B,lspan);
      [~,F1_lam]=Fsph(B,1-lspan);
      Th=(E1*Flam+E2*F1_lam)/E;
      Tspan=Th+Tp+Teff;
    elseif id=='Num'
      dTQ = @(TQ,r)[-TQ(2)/k/Acond(r); Pconv(r)*havg*(Teff-TQ(1))+qgen*Acond(r)];
      Integral=@(Q1)RK4(dTQ,[T(1); Q1],rspan);
      BC=@(TQ) TQ(1,end)-T(2);
      Q1=NSolve(@(Q1)BC(Integral(Q1)),1);
      [TQ, rspan] = Integral(Q1);
      Tspan=TQ(1,:);
    else
      error('plate id must = ''Rec'', ''Bar'', ''Cyl'', ''Sph'' or ''Num''');
    end
  end
end
