function [F12, A1, F21, A2] = disk2cone(rd, rc, h, L)
%[F12, A1, F21, A2] = disk2cone(rd, rc, h, L)
%
%    View factor from a disk to a coaxial cone.
%
%  Inputs:
%
%    rd = radius of disk
%    rc = radius of cone
%    h  = height of cone
%    L  = distance from disk to cone
%
%  Outputs:
%
%    F12 = view factor from disk to cone
%    A1  = area of the disk
%    F21 = view factor from cone to disk
%    A2  = area of the cone
%
%  Reference:
%
%    A Catalog of Radiation Heat Transfer Configuration Factors,
%    third edition, John R. Howell, The University of Texas at Austin
%
%      http://www.thermalradiation.net/sectionc/C-48.html
%
%==========================================================================

alpha = atan(rc/h);  %  Cone angle

A1 = pi*rd^2;
l  = h/cos(alpha);
A2 = pi*rc*l;

S = L/rd;
R = rc/rd;
X = S + R*cot(alpha);
A = sqrt(X^2 + (1 + R)^2);
B = sqrt(X^2 + (1 - R)^2);
C = sqrt(cos(alpha) + S*sin(alpha));
D = sqrt(cos(alpha) - S*sin(alpha));
E = R*cot(alpha) - S;

if alpha > atan(1/S)

  F12 = (R^2 + X^2 + 1 - sqrt((1 + R^2 + X^2)^2 - 4*R^2))/2;
  
else
  
  F12 = (-A*B*atan((A*C)/(B*D)) + (1 + S^2)*atan(C/D) + ...
        (sin(alpha)/cos(alpha)^2)*(X*E*atan((C*D)/X) + ...
        + S^2*atan((D*D/S)) + (C*D)^2*(atan(X/(C*D)) - atan(S/(C*D)))) + ...
        ((R*(X + S))/sin(2*alpha) - S*R*tan(alpha))*acos(-S*tan(alpha)))/pi;

end

F21 = F12*A1/A2;  %  F12*A1 = F21*A2 reciprocity
