function D = dmatr(L,fi,theta,psi)
%
% Function dmatr calculates the D-matrices.
%
% Version 1.0, 060412
%
fid2 = fopen('prodatsres.m', 'a');
fprintf(fid2,'DM-ANGLES  %7.4f %7.4f %7.4f\n', fi, theta, psi);
T = theta/2;
cost = cos(T);
sint = sin(T);
N = 2*L + 1;
for K1 = 1:N
  M1 = K1 - L - 1;
  for K2 = 1:N
    M2 = K2 - L - 1;
    L1 = 1;
    if (M1 - M2) > 0
      L1 = M1 - M2 + 1;
    end
    L2 = L + M1 + 1;
    if (M1 + M2) > 0
      L2 = L - M2 + 1;
    end
    D(K1,K2) = 0;
    for I1 = L1:L2
      I = I1 - 1;
      T = 1;
      K3 = L - M1;
      K4 = L + M1;
      K5 = L - M2;
      K6 = L + M2;
      T = T*fac(K3)*fac(K4)*fac(K5)*fac(K6);
      K3 = L - M2 - I;
      K4 = L + M1 - I;
      K5 = I + M2 - M1;
      T = sqrt(T)/(fac(K3)*fac(K4)*fac(K5)*fac(I));
      N2 = (-1)^I;
      T = T*N2;
      K3 = 2*L - 2*I - M2 + M1;
      Q1 = cost^K3;
      K4 = 2*I + M2 - M1;
      Q2 = sint^K4;
      T = T*Q1*Q2;
      D(K1,K2) = D(K1,K2) + T;
    end
    afi = -i*M1*fi;
    apsi = -i*M2*psi;
    D(K1,K2) = D(K1,K2)*exp(afi+apsi);
  end
end
for K1 = 1:N
  for K2 = 1:N
    Dord(1,K1,K2) = real(D(K1,K2));
    Dord(2,K1,K2) = imag(D(K1,K2));
  end
end
for K1 = 1:N
  fprintf(fid2,'(%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f) (%7.4f %7.4f)\n', Dord(1:2,K1,1:N));
  fprintf(fid2,'\n');
end
fclose(fid2);
return