// Devil Plate by Abderrahman Taha 15/10/2019 (MathMod version 12)
/*

Cst:
c = 1/1000
Th = 1/10

X min/max -+45/10
Y min/max -25/10;+45/10 
Z min/max -+75/10 // Z is up!
*/

#version 3.7;

global_settings { assumed_gamma 2.2 }

#include "colors.inc"
#include "textures.inc"
#include "metals.inc"
#include "functions.inc"


#declare H = 1e-4;
//#declare H2 = 2*H;
#declare Th = 1/10;

#declare Xmin = -45/10;
#declare Xmax = 45/10;
#declare Ymin = -25/10;
#declare Ymax = 45/10;
#declare Zmin = -75/10;
#declare Zmax = 75/10;

// Devil=x^4+2*x^2*z^2-(36/100)*x^2-y^4+(25/100)*y^2+z^4
#declare DevilFn = function {
	pow(x,4)+2*pow(x,2)*pow(z,2)-(36/100)*pow(x,2)-pow(y,4)+(25/100)*pow(y,2)+pow(z,4)
	}

// Devil2=Devil(x,sqrt(y*y+z*z)-(15/10),z,t) 
#declare DevilFn2 = function(x,y,z) {
	DevilFn (
		x,
		sqrt(y*y+z*z)-(15/10),
		z
	       )	
	}

// IsoExterior=Devil2(x,y,sqrt(x*x+z*z)-(15/10),t)
#declare IsoExt = function(x,y,z) {
	DevilFn2 (
		x,
		y,
		sqrt(x*x+z*z)-(15/10)
	       )	
	}

// DFx=((IsoExterior(x+c,y,z,t)-IsoExterior(x,y,z,t))/c)
#declare DFX = function { (IsoExt(x+H,y,z)-IsoExt(x,y,z))/H }

// DFy=((IsoExterior(x,y+c,z,t)-IsoExterior(x,y,z,t))/c)
#declare DFY = function { (IsoExt(x,y+H,z)-IsoExt(x,y,z))/H }

// DFz=((IsoExterior(x,y,z+c,t)-IsoExterior(x,y,z,t))/c)
#declare DFZ = function { (IsoExt(x,y,z+H)-IsoExt(x,y,z))/H }

// R=x/sqrt(x*x+y*y+z*z)
#declare R = function { x/sqrt(x*x+y*y+z*z) }

// Iso=(IsoExterior(x+t*Th*R(DFx(x,y,z,t),DFy(x,y,z,t),DFz(x,y,z,t),t),y+t*Th*R(DFy(x,y,z,t),DFz(x,y,z,t),DFx(x,y,z,t),t),z+t*Th*R(DFz(x,y,z,t),DFx(x,y,z,t),DFy(x,y,z,t),t),t))
#declare Iso = function(x,y,z) {
	IsoExt (
		x+Th*R(DFX(x,y,z),DFY(x,y,z),DFZ(x,y,z)),
		y+Th*R(DFY(x,y,z),DFZ(x,y,z),DFX(x,y,z)),
		z+Th*R(DFZ(x,y,z),DFX(x,y,z),DFY(x,y,z))
		)	
        }

// ThickIsoExterior=(Iso(x,y,z,1)*Iso(x,y,z,-1))
#declare ThickIsoExterior = function(x,y,z) { (Iso(x,y,z)*Iso(x,y,z)) }

// F(x,y,z,t) = if(y<(44/10)&y>-(24/10),ThickIsoExterior((10*x/15),y,z,t),1)
#declare DP_Fn = function(x,y,z) { 
select ( 
        y<(44/10) & y>-(24/10),
        1,
        ThickIsoExterior((10*x/15),y,z)
        
       )
}

//#declare Devil_Plate = 
isosurface {
  function {DP_Fn(x,y,z)

  }
        contained_by{box{<Xmin,Ymin,Zmin>, <Xmax, Ymax, Zmax>}}
        threshold 0.01 // 1.0 Trace Time 1 h 38 min
        accuracy 0.1   // 0.01
        max_gradient 2352446 
        open
        
        pigment {Yellow}
	
	rotate 90*y
	rotate 0*z // -90*z
	rotate 0*x  // 90*x
}   

#declare Win = image_width/image_height;

background { color <0.0,0.0,0.0> }

camera {
    location <0, 10, -15> // faster at location <0, 10, -25>
    up y
    right Win*x
    look_at <0.0, 0.0, 0.0>
    angle 54
}

light_source {
    <7, 25, -35>
    color White
}

/*
Enclosing the inequalities in parentheses evaluates them as a Boolean.
If you multiply them, if either one is false, then the product is zero.
If they are both true, then you get 1 * 1 = 1;
so you could try:

F(x,y,z) =
select ( 
        (y<(44/10)) * (y>-(24/10)),
        1, 
        ThickIsoExterior((10*x/15),y,z)
)

A zero gets added in to use the "4-parameter version" of select, and then since
false = 0, the 1 is the next parameter, and true = 1, so ThickIsoExterior is the
last parameter.
*/