// intersection of two kummer surfaces...

width=600;
height=600;

double  Pi=2*arccos(0);

polyxyz p=z-1+sqrt(2)*x;
polyxyz q=z-1-sqrt(2)*x;
polyxyz r=z+1+sqrt(2)*y;
polyxyz s=z+1-sqrt(2)*y;

double  sf=0.25;
scale_x=sf; scale_y=sf; scale_z=sf;

rot_x=1.7453; rot_y=0.0; rot_z=0.0;

// rot_x=2.4238; rot_y=5.8424; rot_z=4.8279;

double  mu=1.3;
double  lambda=(3*mu-1)/(3-mu);

double  parameter=0.789;

polyxyz sph=x^2+y^2+z^2-mu;

polyxyz contact=1;
polyxyz kummer=sph^2-lambda*p*q*r*s;

kummer = rotate( kummer,
                 Pi/2-arccos(sqrt(2/3)),
                 yAxis);

illumination = ambient_light+
               diffuse_light+
               reflected_light+
               transmitted_light;

light1_x    =100;
light1_y    =100;
light1_z    =100;
light1_red  =255;
light1_green=255;
light1_blue =255;
light1_vol  = 80;

light3_x    =-100;
light3_y    =100;
light3_z    =100;
light3_red  =255;
light3_green=255;
light3_blue =255;
light3_vol  = 60;

light8_x    =  0;
light8_y    =-30;
light8_z    =100;
light8_red  =255;
light8_green=255;
light8_blue =255;
light8_vol  = 40;

transparence=0;
transmitted=0;
smoothness=40;
thickness=0;

transparence2=15;
diffuse2=60;
reflected2=80;
transmitted2=60;
smoothness2=40;
thickness2=11;

surface_red  =40;  inside_red  =surface_red;
surface_green=100;  inside_green=surface_green;
surface_blue =255;  inside_blue =surface_blue;

surface2_red  =255;  inside2_red  =surface2_red;
surface2_green=10;  inside2_green=surface2_green;
surface2_blue =10;  inside2_blue =surface2_blue;

antialiasing=4;
dither_colors=yes;
optimize_colors=no;
normalize_brightness=no;

perspective=central;
spec_z=25.0;
radius=9.0;

int num_images=1500;
int i=0;
double a=10.0;

    a=4.8*((1.0*i)/num_images) + 0.75;

    parameter=(1.0/1907712.0)*
            (-47692800.0+19077120.0*a
            -4226313600.0*a^2+23686040640.0*a^3
            -55926328924.0*a^4+74754566428.0*a^5
            -63384608867.0*a^6+36012423848.0*a^7
            -14056816929.0*a^8+3784317042.0*a^9
            -690880901.0*a^10+81687188.0*a^11
            -5643163.0*a^12+172918.0*a^13)/
                (a^2*(a-5)^2);

    contact=parameter^2*p*q+
        2*parameter*(1-parameter)*sph+
        (1-parameter)^2*lambda*r*s;

    surface=rotate( kummer,
                 i*(2*Pi)/num_images+Pi/4,zAxis );
    contact=rotate( contact,
                 Pi/2-arccos(sqrt(2/3)),
                 yAxis);
    surface2=rotate( contact,
                 i*(2*Pi)/num_images+Pi/4,zAxis );
    draw_surface;

    cutsurface1=surface2;
    cut_with_surface;