// --------------------------------
//  This is a quintic surface with
//  31 double points and hence
//  a Togliatti surface.
// --------------------------------

// useful constants

double  Pi=2*arccos(0);

first  = translate;
second = scale;
third  = rotate;

// fix rotation

rot_x=Pi/2+0.3;
rot_y=0.0;
rot_z=0.0;

// fix translation

origin_z = 0.6;
spec_z = 35.0;
perspective=central;

// fix clipping area

center_y=0.0;
center_z=0.0;
radius=8.0;

// fix scaling

double  sf = 0.35;

scale_x=sf;
scale_y=sf;
scale_z=sf;

// fix illumination

ambient=20;
smoothness=25;
normalize_brightness=yes;
antialiasing=4;

light1_x=-100;
light1_y=100;
light1_z=100;
light1_vol=100;

light2_x=100;
light2_y=20;
light2_z=50;
light2_vol=60;

light3_x=-100;
light3_y=100;
light3_z=100;
light3_vol=0;

surface_red=120;
surface_green=200;
surface_blue=100;
inside_red=surface_red-20;
inside_green=surface_green-20;
inside_blue=surface_blue-20;

// generate surface

double  u1=(sqrt(5)-1)/4;
double  u2=(sqrt(5)+1)/4;
double  v1=(5+sqrt(5))/8;
double  v2=(5-sqrt(5))/8;

double  c=sqrt(5-sqrt(5))/2;
double  a=-(8/5)*(1+1/sqrt(5))*sqrt(5-sqrt(5));
double  r=(1+3*sqrt(5))/4;

poly q=(1-c*z)*(x^2+y^2-1+r*z^2)^2;

poly f=(x-z)*((u1*x-z)^2-v1*y^2)*((u2*x+z)^2-v2*y^2);

poly h1=x-z;
poly h2=cos(2*Pi/5)*x-sin(2*Pi/5)*y-z;
poly h3=cos(4*Pi/5)*x-sin(4*Pi/5)*y-z;
poly h4=cos(6*Pi/5)*x-sin(6*Pi/5)*y-z;
poly h5=cos(8*Pi/5)*x-sin(8*Pi/5)*y-z;

poly F=h1*h2*h3*h4*h5;

poly quintic=a*F+q;

quintic = rotate( quintic,0.5,zAxis );
surface = quintic;

// draw surface

clear_screen;
draw_surface;