You are doing everything right, in principle. The heart of the problem is not related to `tikz-3dplot-circleofsphere` nor `3dtools`. Rather, it is the inaccurracies that come from computations in LaTeX that are to blame. In general, in LaTeX computations can be done by playing with lengths or by using more elaborate tools like the `fpu` library or `xfp`. The results of these methods differ. In the situation at hand, a tiny difference has large impacts since one amplifies them. The radius of the circle and the sphere coincide, at least in theory. You compute two different radii. According to the computation (I switched to `fpu` but it does not solve the problem), the sphere has a radius `R=6.428496000000000` and the circle, i.e. distance between `I` and `S` is `R_S=6.42527000000000`. The relative error is as small as `Delta R=0.00050183`. Nonethless, the `acos` of the ratio of the radii is, according to pgf, as large as `-1.28128`. (The true value is even larger, the trigonometric functions are known to be very inaccurate in pgf when the arguments are very close to the special vaues `0`, `90 degrees` and so on, but this is not really the problem here.) So the upshot is that with `pgf` you will get an incorrect result, at least proceeding this way. You can just set the gamma angle to 0 to obtain the correct result. It would be interesting to know if `xfp` does much better here. 

The problem disappears when one does the computations completely with `fpu`. I updated the `3dtools` library accordingly, the result can be found at https://github.com/marmotghost/tikz-3dtools. Apart from solving this issue, the updated library is arguably also cleaner since it does not introduce all the lengths that were introduced in the previous library. Notice, though that in general the problem persists, see e.g. https://github.com/pgf-tikz/pgf/issues/874 (and the reaction by the maintainer of Ti*k*Z who was informed about this multiple times before the issue got reported). Let me mention that it is in principle possible to solve these problems, along with all the `dimension too large` errors once and for all. 

Note however that the computation of the gamma angle in this way is not too well motivated either. You get the elevation of the circle by comparing radii. However, an arguably much cleaner determination is to compute the distance of the centers of the circle and the sphere, and to compute its `asin` to obtain the gamma angle. In theory, this distance is 0, but even with the slight offset the result is correct. 

```
\documentclass[border=2mm,12pt,tikz]{standalone}
\usepackage{tikz-3dplot-circleofsphere}
\usetikzlibrary{calc}
\usetikzlibrary{fpu}
\def\pgfmathsetmacroFPU#1#2{\begingroup%
	\pgfkeys{/pgf/fpu,/pgf/fpu/output format=fixed}%
	\pgfmathsetmacro{#1}{#2}%
	\pgfmathsmuggle#1\endgroup}%
\usetikzlibrary{3dtools} 
\tikzset{intersection of line trough/.code args={#1 and #2 with plane
		containing #3 and normal #4}{%
		\pgfmathsetmacro{\ltest}{abs(TD("#2o#4")-TD("#1o#4"))}%
		\ifdim\ltest pt<0.01pt            
		\message{Plane and line are parallel!^^J}
		\pgfmathsetmacro{\myd}{0}
		\else
		\pgfmathsetmacro{\myd}{(TD("#3o#4")-TD("#1o#4"))/(TD("#2o#4")-TD("#1o#4"))}%
		\fi
		\pgfmathsetmacro{\myP}{TD("#1+\myd*#2-\myd*#1")}%
		\pgfkeysalso{insert path={%
				(\myP)
		}}
}}

\begin{document}
	
\tdplotsetmaincoords{70}{100}
\begin{tikzpicture}[scale=1,tdplot_main_coords,declare function={a=7;b=7;c=7;h=10;}]
	\path (0,0,0) coordinate (A)
	(c,0,0) coordinate (B)
	[/pgf/fpu,/pgf/fpu/output format=fixed]
	({(pow(b,2) + pow(c,2) - pow(a,2))/(2*c)},{sqrt((a+b-c) *(a-b+c) *(-a+b+c)* (a+b+c))/(2*c)},0) coordinate (C)
	(0,0,h) coordinate (S);
	
	\pic[draw=none,/pgf/fpu,/pgf/fpu/output format=fixed]{3d circle through 3 points={A={(A)},B={(B)},C={(C)},center name=O}};
	\path[overlay] [3d coordinate={(H)=0.5*(S)+0.5*(A)},
	3d coordinate={(n)=(S)-(A)},
	3d coordinate={(D)=2*(O)-(A)},
	3d coordinate={(n1)=(A)-(B)x(A)-(C)},3d coordinate={(T)=(O) + (n1)} ];
	
	\path[intersection of line trough={(O) and (T) with plane containing (H) and normal (n)}] coordinate (I);
	\pgfmathsetmacro{\R}{sqrt(TD("(S)-(I)o(S)-(I)")}
	\pgfmathsetmacroFPU{\Rfpu}{sqrt(TD("(S)-(I)o(S)-(I)")}
	\pgfmathsetmacro{\k}{sqrt(TD("(A)-(O)o(A)-(O)")}
	\begin{scope}[tdplot_screen_coords]
	\draw[blue, thick] (I) circle[radius=\R];
	\end{scope}
	\pic[draw=none,/pgf/fpu,/pgf/fpu/output format=fixed]{3d circle through 3 points={A={(A)},B={(D)},C={(S)},center name=T'}};
	
	\begin{scope}[shift={(I)}]
	\path[overlay] [3d coordinate={(myn)=(S)-(D)x(S)-(A)},
	3d coordinate={(S-T')=(S)-(T')},
	3d coordinate={(I-T')=(I)-(T')}];
	\pgfmathsetmacro{\myaxisangles}{axisangles("(myn)")}
	\pgfmathsetmacro{\myalpha}{{\myaxisangles}[0]}
	\pgfmathsetmacro{\mybeta}{{\myaxisangles}[1]}
	\pgfmathsetmacroFPU{\RSfpu}{sqrt(TD("(S-T')o(S-T')")}
	\pgfmathsetmacroFPU{\Rerror}{(\Rfpu-\RSfpu)/\Rfpu}
	\pgfmathsetmacroFPU{\mygamma}{-acos(\RSfpu/\Rfpu)}
	\typeout{R_S=c,Delta R=\Rerror,fake gamma=\mygamma}
	%
	\pgfmathsetmacro{\mydist}{sqrt(TD("(I-T')o(I-T')"))}
	\pgfmathsetmacroFPU{\mygamma}{asin(\mydist/\Rfpu)}
	\typeout{\mydist,\mygamma}
	\tdplotCsDrawCircle[tdplotCsFront/.style={thick,purple}]{\R}{\myalpha}{\mybeta}{\mygamma} 
	\end{scope}
	\foreach \p in {A,B,C,S,O,I,D,T'}
	\draw[fill=black] (\p) circle (1.5 pt);
	\foreach \p/\g in {A/-90,B/-90,C/-90,S/150,O/90,I/0,D/90}
	\path (\p)+(\g:3mm) node{$\p$};
	\draw[dashed] (S) -- (A)  (S) -- (B)  (S) -- (C) (A) -- (B) -- (C)  --cycle (S) -- (D) (A) -- (D);
	\end{tikzpicture}
\end{document}
```

![Screen Shot 2020-08-14 at 11.23.28 AM.png](/image?hash=f40890b7e870efba9b82c4c35288725788a3cb21891cd0e8a6fedd84fcd42a71)