Surf plot with tikz

Question

I am trying to regenerate a surf plot from Matlab using 3-d-plot of tikz.

The graph should look like this:

But with my code:

\documentclass[]{article}

\usepackage{pgfplots} 
\usepackage{pgfgantt}
\usepackage{pdflscape}
\pgfplotsset{compat=newest} 
\pgfplotsset{plot coordinates/math parser=false}

\usepackage{tikz}
\usepackage{tikz-3dplot}
\usetikzlibrary{intersections,fadings,decorations.pathreplacing}
 \usepgfplotslibrary{patchplots}

\begin{document}

\begin{tikzpicture}
\begin{axis}[grid=major,
colormap/viridis]
\addplot3[surf, shader=faceted, samples =80, domain=0:6.28, y domain = 0:6.28]
{(45069.0*(0.000019484*sin(y) + 0.000008086700*sin(y)*sin(y)*sin(y)- 0.0000012105*x*cos(y)))/(0.0031581*cos(2.0*y) - 0.007735)};
\end{axis}
\end{tikzpicture}
\end{document}

I am getting this output:

I can not figure out, what the problem is. Can you help me? Find attached also the matlab live script generating the required output:

x = sym('x_%d', [4, 1], {'real'}); % system states
u = sym('u_%d', [1, 1], {'real'}); % control input

x_1_val_vec = 0:0.1:2*pi;
x_2_val_vec = 0:0.1:2*pi;
[x_grid, y_grid] = meshgrid(x_1_val_vec, x_2_val_vec);
%f_subs_plot = subs(f_subs(4), {x(3), x(4), u(1)}, {0, 0, 0});
f_subs_plot= (45069.0*(0.000019484*sin(x(2)) + 8.0867e-6*sin(x(2))^3 - 1.2105e-6*x(1)*cos(x(2))))/(0.0031581*cos(2.0*x(2)) - 0.007735);

f_con_eval = eval(subs(f_subs_plot, {x(1), x(2)}, {x_grid, y_grid}));

figure()
h_fig = surf(x_grid, y_grid, f_con_eval, 'LineWidth', 0.1);
h_fig_ax = get(h_fig, 'Parent');
hold on
z_zero = zeros(size(x_grid, 1));
zero_fig = surf(x_grid, y_grid, z_zero, 'FaceAlpha',0.5, 'LineStyle', 'none');

xlabel('x_1')
ylabel('x_2')
zlabel('g')

Answer 1

The main part of this question is a duplicate of e.g. this one, and the answer is: "use trig format=rad". However, there is an interesting side aspect: how can one add the plane. This can be done by letting TikZ find the intersection points of the outer boundaries of the surface numerically (in a \savebox) and then drawing the plot step by step to have 3d ordering.

\documentclass{article}
\usepackage{pgfplots} 
\pgfplotsset{compat=1.16} 
\usetikzlibrary{intersections,calc}
\newsavebox\Schummel
\begin{document}
\begingroup\tikzset{declare function={f(\x,\y)=(45069.0*(0.000019484*sin(\y) +
0.000008086700*sin(\y)*sin(\y)*sin(\y)-
0.0000012105*\x*cos(\y)))/(0.0031581*cos(2.0*\y) - 0.007735);}}
\savebox\Schummel{\begin{tikzpicture}
 \begin{axis}[trig format=rad]
 \addplot[name path=plot1,domain=0:{2*pi}] {f(0,x)};
 \addplot[name path=plot2,domain=0:{2*pi}] {f(2*pi,x)};
 \draw[name path global=axis] (-0.1,0) -- (2*pi+0.1,0);
 \path (0,0) coordinate (O) (1,0) coordinate (X);
 \path[name intersections={of=axis and plot2,total=\t}]
  foreach \X in {1,...,\t} {let \p1=(intersection-\X),\p2=($(X)-(O)$) in
  \pgfextra{\pgfmathsetmacro{\myx}{\x1/\x2}
  \ifnum\X=1   
   \xdef\lstX{\myx}
  \else
   \xdef\lstX{\lstX,\myx}
  \fi}};
 \end{axis} 
\end{tikzpicture}}
\begin{tikzpicture}
\begin{axis}[grid=major,trig format=rad,view={20}{60},colormap/viridis,
zlabel={$f_4$},zlabel style={rotate=-90},
samples=81,domain = 0:{2*pi}]
 \addplot3[surf, shader=faceted,samples y=21,domain y=0:pi/2] {f(x,y)};
 \pgfmathsetmacro{\yone}{{\lstX}[0]}
 \pgfmathsetmacro{\ytwo}{{\lstX}[1]}
 \path[opacity=0.2,fill=cyan]  (0,0,0) -- (2*pi,\yone,0) -- (2*pi,\ytwo,0)
 -- (0,pi,0) -- cycle;
 \addplot3[surf, shader=faceted,samples y=61,domain y=pi/2:2*pi] {f(x,y)};
 \path[opacity=0.2,fill=cyan]   (0,pi,0) -- (2*pi,\ytwo,0)
 -- (2*pi,2*pi,0) -- (0,2*pi,0) -- cycle;
\end{axis}
\end{tikzpicture}
\endgroup
\end{document}

If you change the view angles, you may have to add the pieces in a different order.