END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

double Icon(double x, double z, double y){
double rde;
const int Maxiter=8;
double MaxRad=60.0;
double p2,p3,p4,p5,p6;
double rad=sqrt(x*x+y*y+z*z);
int op;

op=4;// Four Sides
double Alpha=-1.327888;
double Beta=-1.116979;
double Gamma=1.364337;
double Lambda=1.722694;

int op1=op-2;
int op2=op-1;
double p3z;
int Counter=0;
double xr,yr,zr,nx,ny,nz;
do{
xr=x;
yr=y;
zr=z;
const double p1rad=sqrt(x*x+y*y+z*z); //the standard symmetric formula doesn't use the square root, so you can remove it too
const double angle=atan2(y,x)*op2;
const double angle1=atan2(sqrt(x*x+y*y),z)*op2;
const double radius=pow(p1rad,op2);
const double sinp=sin(angle1)*radius;
p2=cos(angle)*sinp;
p3=sin(angle)*sinp;
p3z=cos(angle1)*radius;
p5=x*p2-y*p3-z*p3z;
p6=Lambda+Alpha*p1rad+Beta*p5;
x=p6*x+Gamma*p2;
y=p6*y-Gamma*p3;
z=p6*z-Gamma*p3z;
Counter++;
rad=sqrt(x*x+y*y+z*z);
if(rad>MaxRad)break;
}while(Counter<Maxiter);
const double minr=sqrt(xr*xr+yr*yr+zr*zr);
const double potential=sqrt(minr/MaxRad);
double itential=(double(Maxiter-Counter)/double(Maxiter));
double itential1=(double(Maxiter-Counter+1)/double(Maxiter));
const double potential1=1-sqrt(MaxRad/rad);
const double totalpot=(potential*0.5+potential1*0.5);
double potiter=(itential1*totalpot+itential*(1-totalpot));
 rde=pow(potiter,Counter);
return rde;
}