To draw the circle, we need to switch to the plane in which B, C and E (and here also F) sit, and to compute the center and the radius of the circle.
This answer focuses on how to transform the coordinate system (since you seem to know the center and the radius and the question how to get that coordinate system has a surprisingly simple answer).
Call the normalized normal vector of the plane n. We need to find rotation angles in such a way that the rotated z axis coincides with n. However, the rotated z axis is simply
D.(0,0,1)
where the rotation matrix D is given on p. 7 of the tikz-3dplot manual
So we need to solve
(0) (nx)
D |0| = |ny|
(1) (nz)
which has the solution
beta = arccos(nz)
alpha = arccos(nx/sin(beta))
For the users convenience, the circles can be inserted via a style
\draw[red,thick,circle in plane with normal={{\mynormal} with radius {\r} around (I)}];
This style also comes with a correction. The earlier version of the answer sometimes yielded wrong output, with the reason being as simple as a sign reflecting that cos(x)=cos(-x) (GRRRR). Neither the centers nor the radii of the circles are computed here, I got them all out of our chat. Sorry for the sign mistake!
\documentclass[tikz,border=3.14mm]{standalone}
\usepackage{tikz-3dplot}
\usetikzlibrary{3d}
\usetikzlibrary{calc}
\makeatletter
\newcounter{smuggle}
\DeclareRobustCommand\smuggleone[1]{%
\stepcounter{smuggle}%
\expandafter\global\expandafter\let\csname smuggle@\arabic{smuggle}\endcsname#1%
\aftergroup\let\aftergroup#1\expandafter\aftergroup\csname smuggle@\arabic{smuggle}\endcsname
}
\DeclareRobustCommand\smuggle[2][1]{%
\smuggleone{#2}%
\ifnum#1>1
\aftergroup\smuggle\aftergroup[\expandafter\aftergroup\the\numexpr#1-1\aftergroup]\aftergroup#2%
\fi
}
\makeatother
\def\parsecoord(#1,#2,#3)>(#4,#5,#6){%
\def#4{#1}%
\def#5{#2}%
\def#6{#3}%
\smuggle{#4}%
\smuggle{#5}%
\smuggle{#6}%
}
\def\SPTD(#1,#2,#3).(#4,#5,#6){((#1)*(#4)+1*(#2)*(#5)+1*(#3)*(#6))}
\def\VPTD(#1,#2,#3)x(#4,#5,#6){((#2)*(#6)-1*(#3)*(#5),(#3)*(#4)-1*(#1)*(#6),(#1)*(#5)-1*(#2)*(#4))}
\def\VecMinus(#1,#2,#3)-(#4,#5,#6){(#1-1*(#4),#2-1*(#5),#3-1*(#6))}
\def\VecAdd(#1,#2,#3)+(#4,#5,#6){(#1+1*(#4),#2+1*(#5),#3+1*(#6))}
\newcommand{\RotationAnglesForPlaneWithNormal}[5]{%\typeout{N=(#1,#2,#3)}
\foreach \XS in {1,-1}
{\foreach \YS in {1,-1}
{\pgfmathsetmacro{\mybeta}{\XS*acos(#3)}
\pgfmathsetmacro{\myalpha}{\YS*acos(#1/sin(\mybeta))}
\pgfmathsetmacro{\ntest}{abs(cos(\myalpha)*sin(\mybeta)-#1)%
+abs(sin(\myalpha)*sin(\mybeta)-#2)+abs(cos(\mybeta)-#3)}
\ifdim\ntest pt<0.1pt
\xdef#4{\myalpha}
\xdef#5{\mybeta}
\fi
}}
}
\tikzset{circle in plane with normal/.style args={#1 with radius #2 around #3}{
/utils/exec={\edef\temp{\noexpand\parsecoord#1>(\noexpand\myNx,\noexpand\myNy,\noexpand\myNz)}
\temp
\pgfmathsetmacro{\myNx}{\myNx}
\pgfmathsetmacro{\myNy}{\myNy}
\pgfmathsetmacro{\myNz}{\myNz}
\pgfmathsetmacro{\myNormalization}{sqrt(pow(\myNx,2)+pow(\myNy,2)+pow(\myNz,2))}
\pgfmathsetmacro{\myNx}{\myNx/\myNormalization}
\pgfmathsetmacro{\myNy}{\myNy/\myNormalization}
\pgfmathsetmacro{\myNz}{\myNz/\myNormalization}
% compute the rotation angles that transform us in the corresponding plabe
\RotationAnglesForPlaneWithNormal{\myNx}{\myNy}{\myNz}{\tmpalpha}{\tmpbeta}
%\typeout{N=(\myNx,\myNy,\myNz),alpha=\tmpalpha,beta=\tmpbeta,r=#2,#3}
\tdplotsetrotatedcoords{\tmpalpha}{\tmpbeta}{0}},
insert path={[tdplot_rotated_coords,canvas is xy plane at z=0,transform shape]
#3 circle[radius=#2]}
}}
\begin{document}
\foreach \X in {5,15,...,355} % {50}%
{\tdplotsetmaincoords{70}{\X}
\begin{tikzpicture}[tdplot_main_coords,scale=1]
\path [tdplot_screen_coords,use as bounding box] (-7,-3) rectangle (7,7);
\pgfmathsetmacro\a{1}
\pgfmathsetmacro\h{7}
\pgfmathsetmacro\rprime{5*sqrt(\a^2*\h^2/(25*\a^2+\h^2))*(1/2))}
\pgfmathsetmacro\r{(1/2)*sqrt((400*\a^4+9*\a^2*\h^2)/(16*\a^2+\h^2))}
% definitions
\path
coordinate(A) at (0,0,0)
coordinate (B) at (4*\a,0,0)
coordinate (C) at (4*\a,3*\a,0)
coordinate (S) at (0,0,\h)
coordinate (E) at ({4*\h^2*\a/(16*\a^2+\h^2)}, 0, {16*\h*\a^2/(16*\a^2+\h^2)})
coordinate (F) at ({4*\h^2*\a/(25*\a^2+\h^2)}, {3*\h^2*\a/(25*\a^2+\h^2)},{25*\h*\a^2/(25*\a^2+\h^2)})
coordinate (I') at ($(F)!0.5!(A) $)
coordinate (I) at ($(C)!0.5!(E) $);
\begin{scope}
\draw[dashed,thick]
(A) -- (C) ;
\draw[thick]
(S) -- (B) (S)-- (A) -- (B)-- (C) -- cycle (A) --(E) (A) --(F);
\end{scope}
\foreach \point/\position in {A/left,B/below,C/right,S/above,E/right,F/below}
{
\fill[black] (\point) circle (1.5pt);
\node[\position=1.5pt] at (\point) {$\point$};
}
% % store the coordinates of E, A and F in marcros
\parsecoord({4*\h^2*\a/(16*\a^2+\h^2)},0,{16*\h*\a^2/(16*\a^2+\h^2)})>(\myEx,\myEy,\myEz)
\parsecoord(0,0,0)>(\myAx,\myAy,\myAz)
\parsecoord({4*\h^2*\a/(25*\a^2+\h^2)},{3*\h^2*\a/(25*\a^2+\h^2)},{25*\h*\a^2/(25*\a^2+\h^2)})>(\myFx,\myFy,\myFz)
\parsecoord(4*\a,0,0)>(\myBx,\myBy,\myBz)
\parsecoord(4*\a,3*\a,0)>(\myCx,\myCy,\myCz)
% % compute the normal of the plane in which E, B and C sit
\def\mynormal{\VPTD({\myEx-\myAx},{\myEy-\myAy},{\myEz-\myAz})x({\myFx-\myAx},{\myFy-\myAy},{\myFz-\myAz})}
\edef\temp{\noexpand\parsecoord\mynormal>(\noexpand\myNx,\noexpand\myNy,\noexpand\myNz)}
\draw[blue,thick,circle in plane with normal={{\mynormal} with radius {\rprime}
around (I')}];
\def\mynormal{\VPTD({\myEx-\myBx},{\myEy-\myBy},{\myEz-\myBz})x({\myCx-\myBx},{\myCy-\myBy},{\myCz-\myBz})}
\draw[red,thick,circle in plane with normal={{\mynormal} with radius {\r}
around (I)}];
\node[anchor=north east] at (current bounding box.north east) {$\theta=\tdplotmaintheta^\circ,\phi=\tdplotmainphi^\circ$};
\end{tikzpicture}}
\end{document}
Notes: (i) to myself this trick may be important for the analytic discrimination of hidden vs. visible parts of objects. If I rediscover this one day, I can at least claim I knew that before. (ii) whoever tries to further develop this very nice answer may find this useful.
