rendered paste body/******************************************************************************* * fnintern.cpp * * This module implements the the function type used by iso surfaces and * the function pattern. * * This module is based on code by D. Skarda, T. Bily and R. Suzuki. * It includes functions based on code first introduced by many other * contributors. * * from Persistence of Vision Ray Tracer ('POV-Ray') version 3.7. * Copyright 1991-2003 Persistence of Vision Team * Copyright 2003-2009 Persistence of Vision Raytracer Pty. Ltd. * --------------------------------------------------------------------------- * NOTICE: This source code file is provided so that users may experiment * with enhancements to POV-Ray and to port the software to platforms other * than those supported by the POV-Ray developers. There are strict rules * regarding how you are permitted to use this file. These rules are contained * in the distribution and derivative versions licenses which should have been * provided with this file. * * These licences may be found online, linked from the end-user license * agreement that is located at http://www.povray.org/povlegal.html * --------------------------------------------------------------------------- * POV-Ray is based on the popular DKB raytracer version 2.12. * DKBTrace was originally written by David K. Buck. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins. * --------------------------------------------------------------------------- * $File: //depot/povray/smp/source/backend/vm/fnintern.cpp $ * $Revision: #16 $ * $Change: 4715 $ * $DateTime: 2009/02/21 10:00:39 $ * $Author: chrisc $ *******************************************************************************//********************************************************************************* * NOTICE * * This file is part of a BETA-TEST version of POV-Ray version 3.7. It is not * final code. Use of this source file is governed by both the standard POV-Ray * licences referred to in the copyright header block above this notice, and the * following additional restrictions numbered 1 through 4 below: * * 1. This source file may not be re-distributed without the written permission * of Persistence of Vision Raytracer Pty. Ltd. * * 2. This notice may not be altered or removed. * * 3. Binaries generated from this source file by individuals for their own * personal use may not be re-distributed without the written permission * of Persistence of Vision Raytracer Pty. Ltd. Such personal-use binaries * are not required to have a timeout, and thus permission is granted in * these circumstances only to disable the timeout code contained within * the beta software. * * 4. Binaries generated from this source file for use within an organizational * unit (such as, but not limited to, a company or university) may not be * distributed beyond the local organizational unit in which they were made, * unless written permission is obtained from Persistence of Vision Raytracer * Pty. Ltd. Additionally, the timeout code implemented within the beta may * not be disabled or otherwise bypassed in any manner. * * The following text is not part of the above conditions and is provided for * informational purposes only. * * The purpose of the no-redistribution clause is to attempt to keep the * circulating copies of the beta source fresh. The only authorized distribution * point for the source code is the POV-Ray website and Perforce server, where * the code will be kept up to date with recent fixes. Additionally the beta * timeout code mentioned above has been a standard part of POV-Ray betas since * version 1.0, and is intended to reduce bug reports from old betas as well as * keep any circulating beta binaries relatively fresh. * * All said, however, the POV-Ray developers are open to any reasonable request * for variations to the above conditions and will consider them on a case-by-case * basis. * * Additionally, the developers request your co-operation in fixing bugs and * generally improving the program. If submitting a bug-fix, please ensure that * you quote the revision number of the file shown above in the copyright header * (see the '$Revision:' field). This ensures that it is possible to determine * what specific copy of the file you are working with. The developers also would * like to make it known that until POV-Ray 3.7 is out of beta, they would prefer * to emphasize the provision of bug fixes over the addition of new features. * * Persons wishing to enhance this source are requested to take the above into * account. It is also strongly suggested that such enhancements are started with * a recent copy of the source. * * The source code page (see http://www.povray.org/beta/source/) sets out the * conditions under which the developers are willing to accept contributions back * into the primary source tree. Please refer to those conditions prior to making * any changes to this source, if you wish to submit those changes for inclusion * with POV-Ray. * *********************************************************************************/// frame.h must always be the first POV file included (pulls in platform config)#include "backend/frame.h"#include "backend/vm/fncode.h"#include "backend/vm/fnpovfpu.h"#include "backend/vm/fnintern.h"#include "backend/texture/texture.h"#include "backend/texture/pigment.h"#include "backend/math/vector.h"#include "backend/math/mathutil.h"#include "backend/math/splines.h"#include <algorithm>// this must be the last file included#include "base/povdebug.h"namespace pov{using namespace std;/****************************************************************************** Local preprocessor defines******************************************************************************//* * IMPORTANT NOTICE * Never ever use any PARAM_xxx macro after a call to POVFPU_SetLocal * because POVFPU_SetLocal may change the stack base address and thus * the "ptr" passed to the functions may be invalid after the call! [trf] */#define PARAM_X (ptr[0])#define PARAM_Y (ptr[1])#define PARAM_Z (ptr[2])#define PARAM(index) (ptr[index + 3])#define PARAM_N_X(offset) (ptr[offset])#define PARAM_N_Y(offset) (ptr[offset + 1])#define PARAM_N_Z(offset) (ptr[offset + 2])#define PARAM_N(index,offset) (ptr[index + offset + 3])#define ROT2D(p,d,ang) if (p>0) {x2=sqrt(x2+PARAM_Z*PARAM_Z)- d; th=ang*M_PI_180; \ if (th!=0){ PARAM_X= x2*cos(th)-PARAM_Y*sin(th); PARAM_Y= x2*sin(th)+PARAM_Y*cos(th);} else PARAM_X=x2; \ x2=PARAM_X*PARAM_X; y2=PARAM_Y*PARAM_Y;}/****************************************************************************** Global functions******************************************************************************/DBL f_algbr_cyl1(FPUContext *ctx, DBL *ptr, unsigned int fn); // 0DBL f_algbr_cyl2(FPUContext *ctx, DBL *ptr, unsigned int fn); // 1DBL f_algbr_cyl3(FPUContext *ctx, DBL *ptr, unsigned int fn); // 2DBL f_algbr_cyl4(FPUContext *ctx, DBL *ptr, unsigned int fn); // 3DBL f_bicorn(FPUContext *ctx, DBL *ptr, unsigned int fn); // 4DBL f_bifolia(FPUContext *ctx, DBL *ptr, unsigned int fn); // 5DBL f_blob(FPUContext *ctx, DBL *ptr, unsigned int fn); // 6DBL f_blob2(FPUContext *ctx, DBL *ptr, unsigned int fn); // 7DBL f_boy_surface(FPUContext *ctx, DBL *ptr, unsigned int fn); // 8DBL f_comma(FPUContext *ctx, DBL *ptr, unsigned int fn); // 9DBL f_cross_ellipsoids(FPUContext *ctx, DBL *ptr, unsigned int fn); // 10DBL f_crossed_trough(FPUContext *ctx, DBL *ptr, unsigned int fn); // 11DBL f_cubic_saddle(FPUContext *ctx, DBL *ptr, unsigned int fn); // 12DBL f_cushion(FPUContext *ctx, DBL *ptr, unsigned int fn); // 13DBL f_devils_curve(FPUContext *ctx, DBL *ptr, unsigned int fn); // 14DBL f_devils_curve_2d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 15DBL f_dupin_cyclid(FPUContext *ctx, DBL *ptr, unsigned int fn); // 16DBL f_ellipsoid(FPUContext *ctx, DBL *ptr, unsigned int fn); // 17DBL f_enneper(FPUContext *ctx, DBL *ptr, unsigned int fn); // 18DBL f_flange_cover(FPUContext *ctx, DBL *ptr, unsigned int fn); // 19DBL f_folium_surface(FPUContext *ctx, DBL *ptr, unsigned int fn); // 20DBL f_folium_surface_2d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 21DBL f_glob(FPUContext *ctx, DBL *ptr, unsigned int fn); // 22DBL f_heart(FPUContext *ctx, DBL *ptr, unsigned int fn); // 23DBL f_helical_torus(FPUContext *ctx, DBL *ptr, unsigned int fn); // 24DBL f_helix1(FPUContext *ctx, DBL *ptr, unsigned int fn); // 25DBL f_helix2(FPUContext *ctx, DBL *ptr, unsigned int fn); // 26DBL f_hex_x(FPUContext *ctx, DBL *ptr, unsigned int fn); // 27DBL f_hex_y(FPUContext *ctx, DBL *ptr, unsigned int fn); // 28DBL f_hetero_mf(FPUContext *ctx, DBL *ptr, unsigned int fn); // 29DBL f_hunt_surface(FPUContext *ctx, DBL *ptr, unsigned int fn); // 30DBL f_hyperbolic_torus(FPUContext *ctx, DBL *ptr, unsigned int fn); // 31DBL f_isect_ellipsoids(FPUContext *ctx, DBL *ptr, unsigned int fn); // 32DBL f_kampyle_of_eudoxus(FPUContext *ctx, DBL *ptr, unsigned int fn); // 33DBL f_kampyle_of_eudoxus_2d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 34DBL f_klein_bottle(FPUContext *ctx, DBL *ptr, unsigned int fn); // 35DBL f_kummer_surface_v1(FPUContext *ctx, DBL *ptr, unsigned int fn); // 36DBL f_kummer_surface_v2(FPUContext *ctx, DBL *ptr, unsigned int fn); // 37DBL f_lemniscate_of_gerono(FPUContext *ctx, DBL *ptr, unsigned int fn); // 38DBL f_lemniscate_of_gerono_2d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 39DBL f_mesh1(FPUContext *ctx, DBL *ptr, unsigned int fn); // 40DBL f_mitre(FPUContext *ctx, DBL *ptr, unsigned int fn); // 41DBL f_nodal_cubic(FPUContext *ctx, DBL *ptr, unsigned int fn); // 42DBL f_odd(FPUContext *ctx, DBL *ptr, unsigned int fn); // 43DBL f_ovals_of_cassini(FPUContext *ctx, DBL *ptr, unsigned int fn); // 44DBL f_paraboloid(FPUContext *ctx, DBL *ptr, unsigned int fn); // 45DBL f_parabolic_torus(FPUContext *ctx, DBL *ptr, unsigned int fn); // 46DBL f_ph(FPUContext *ctx, DBL *ptr, unsigned int fn); // 47DBL f_pillow(FPUContext *ctx, DBL *ptr, unsigned int fn); // 48DBL f_piriform(FPUContext *ctx, DBL *ptr, unsigned int fn); // 49DBL f_piriform_2d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 50DBL f_poly4(FPUContext *ctx, DBL *ptr, unsigned int fn); // 51DBL f_polytubes(FPUContext *ctx, DBL *ptr, unsigned int fn); // 52DBL f_quantum(FPUContext *ctx, DBL *ptr, unsigned int fn); // 53DBL f_quartic_paraboloid(FPUContext *ctx, DBL *ptr, unsigned int fn); // 54DBL f_quartic_saddle(FPUContext *ctx, DBL *ptr, unsigned int fn); // 55DBL f_quartic_cylinder(FPUContext *ctx, DBL *ptr, unsigned int fn); // 56DBL f_r(FPUContext *ctx, DBL *ptr, unsigned int fn); // 57DBL f_ridge(FPUContext *ctx, DBL *ptr, unsigned int fn); // 58DBL f_ridged_mf(FPUContext *ctx, DBL *ptr, unsigned int fn); // 59DBL f_rounded_box(FPUContext *ctx, DBL *ptr, unsigned int fn); // 60DBL f_sphere(FPUContext *ctx, DBL *ptr, unsigned int fn); // 61DBL f_spikes(FPUContext *ctx, DBL *ptr, unsigned int fn); // 62DBL f_spikes_2d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 63DBL f_spiral(FPUContext *ctx, DBL *ptr, unsigned int fn); // 64DBL f_steiners_roman(FPUContext *ctx, DBL *ptr, unsigned int fn); // 65DBL f_strophoid(FPUContext *ctx, DBL *ptr, unsigned int fn); // 66DBL f_strophoid_2d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 67DBL f_superellipsoid(FPUContext *ctx, DBL *ptr, unsigned int fn); // 68DBL f_th(FPUContext *ctx, DBL *ptr, unsigned int fn); // 69DBL f_torus(FPUContext *ctx, DBL *ptr, unsigned int fn); // 70DBL f_torus2(FPUContext *ctx, DBL *ptr, unsigned int fn); // 71DBL f_torus_gumdrop(FPUContext *ctx, DBL *ptr, unsigned int fn); // 72DBL f_umbrella(FPUContext *ctx, DBL *ptr, unsigned int fn); // 73DBL f_witch_of_agnesi(FPUContext *ctx, DBL *ptr, unsigned int fn); // 74DBL f_witch_of_agnesi_2d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 75DBL f_noise3d(FPUContext *ctx, DBL *ptr, unsigned int fn); // 76DBL f_pattern(FPUContext *ctx, DBL *ptr, unsigned int fn); // 77DBL f_noise_generator(FPUContext *ctx, DBL *ptr, unsigned int fn); // 78DBL f_mandelbulb(FPUContext *ctx, DBL *ptr, unsigned int fn); // 79void f_pigment(FPUContext *ctx, DBL *ptr, unsigned int fn, unsigned int sp); // 0void f_transform(FPUContext *ctx, DBL *ptr, unsigned int fn, unsigned int sp); // 1void f_spline(FPUContext *ctx, DBL *ptr, unsigned int fn, unsigned int sp); // 2/****************************************************************************** Global variables******************************************************************************/const Trap POVFPU_TrapTable[] ={ { f_algbr_cyl1, 5 + 3 }, // 0 { f_algbr_cyl2, 5 + 3 }, // 1 { f_algbr_cyl3, 5 + 3 }, // 2 { f_algbr_cyl4, 5 + 3 }, // 3 { f_bicorn, 2 + 3 }, // 4 { f_bifolia, 2 + 3 }, // 5 { f_blob, 5 + 3 }, // 6 { f_blob2, 4 + 3 }, // 7 { f_boy_surface, 2 + 3 }, // 8 { f_comma, 1 + 3 }, // 9 { f_cross_ellipsoids, 4 + 3 }, // 10 { f_crossed_trough, 1 + 3 }, // 11 { f_cubic_saddle, 1 + 3 }, // 12 { f_cushion, 1 + 3 }, // 13 { f_devils_curve, 1 + 3 }, // 14 { f_devils_curve_2d, 6 + 3 }, // 15 { f_dupin_cyclid, 6 + 3 }, // 16 { f_ellipsoid, 3 + 3 }, // 17 { f_enneper, 1 + 3 }, // 18 { f_flange_cover, 4 + 3 }, // 19 { f_folium_surface, 3 + 3 }, // 20 { f_folium_surface_2d, 6 + 3 }, // 21 { f_glob, 1 + 3 }, // 22 { f_heart, 1 + 3 }, // 23 { f_helical_torus, 10 + 3 }, // 24 { f_helix1, 7 + 3 }, // 25 { f_helix2, 7 + 3 }, // 26 { f_hex_x, 1 + 3 }, // 27 { f_hex_y, 1 + 3 }, // 28 { f_hetero_mf, 6 + 3 }, // 29 { f_hunt_surface, 1 + 3 }, // 30 { f_hyperbolic_torus, 3 + 3 }, // 31 { f_isect_ellipsoids, 4 + 3 }, // 32 { f_kampyle_of_eudoxus, 3 + 3 }, // 33 { f_kampyle_of_eudoxus_2d, 6 + 3 }, // 34 { f_klein_bottle, 1 + 3 }, // 35 { f_kummer_surface_v1, 1 + 3 }, // 36 { f_kummer_surface_v2, 4 + 3 }, // 37 { f_lemniscate_of_gerono, 1 + 3 }, // 38 { f_lemniscate_of_gerono_2d, 6 + 3 }, // 39 { f_mesh1, 5 + 3 }, // 40 { f_mitre, 1 + 3 }, // 41 { f_nodal_cubic, 1 + 3 }, // 42 { f_odd, 1 + 3 }, // 43 { f_ovals_of_cassini, 4 + 3 }, // 44 { f_paraboloid, 1 + 3 }, // 45 { f_parabolic_torus, 3 + 3 }, // 46 { f_ph, 0 + 3 }, // 47 { f_pillow, 1 + 3 }, // 48 { f_piriform, 1 + 3 }, // 49 { f_piriform_2d, 7 + 3 }, // 50 { f_poly4, 5 + 3 }, // 51 { f_polytubes, 6 + 3 }, // 52 { f_quantum, 1 + 3 }, // 53 { f_quartic_paraboloid, 1 + 3 }, // 54 { f_quartic_saddle, 1 + 3 }, // 55 { f_quartic_cylinder, 3 + 3 }, // 56 { f_r, 0 + 3 }, // 57 { f_ridge, 6 + 3 }, // 58 { f_ridged_mf, 6 + 3 }, // 59 { f_rounded_box, 4 + 3 }, // 60 { f_sphere, 1 + 3 }, // 61 { f_spikes, 5 + 3 }, // 62 { f_spikes_2d, 4 + 3 }, // 63 { f_spiral, 6 + 3 }, // 64 { f_steiners_roman, 1 + 3 }, // 65 { f_strophoid, 4 + 3 }, // 66 { f_strophoid_2d, 7 + 3 }, // 67 { f_superellipsoid, 2 + 3 }, // 68 { f_th, 0 + 3 }, // 69 { f_torus, 2 + 3 }, // 70 { f_torus2, 3 + 3 }, // 71 { f_torus_gumdrop, 1 + 3 }, // 72 { f_umbrella, 1 + 3 }, // 73 { f_witch_of_agnesi, 2 + 3 }, // 74 { f_witch_of_agnesi_2d, 6 + 3 }, // 75 { f_noise3d, 0 + 3 }, // 76 { f_pattern, 0 + 3 }, // 77 { f_noise_generator, 1 + 3 }, // 78 { f_mandelbulb, 2 + 3 }, // 79 { NULL, 0 }};const TrapS POVFPU_TrapSTable[] ={ { f_pigment, 0 + 3 }, // 0 { f_transform, 0 + 3 }, // 1 { f_spline, 0 + 1 }, // 2 { NULL, 0 }};const unsigned int POVFPU_TrapTableSize = 80;const unsigned int POVFPU_TrapSTableSize = 3;/****************************************************************************** Functions******************************************************************************/DBL f_algbr_cyl1(FPUContext *ctx, DBL *ptr, unsigned int) // 0{ DBL r, x2, y2, th; x2 = PARAM_X * PARAM_X; y2 = PARAM_Y * PARAM_Y; ROT2D(PARAM(2),PARAM(3),PARAM(4)) PARAM_X = fabs(PARAM_X); r = -(x2 * PARAM_X - x2 + y2); return (-min(PARAM(1), max(PARAM(0) * r, -PARAM(1))));}DBL f_algbr_cyl2(FPUContext *ctx, DBL *ptr, unsigned int) // 1{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, th; ROT2D(PARAM(2),PARAM(3),PARAM(4)) r=-( 2*x2*x2 -3*x2*PARAM_Y +y2 -2*y2*PARAM_Y +y2*y2 ); return( min(PARAM(1), max(PARAM(0)*r,-PARAM(1))) );}DBL f_algbr_cyl3(FPUContext *ctx, DBL *ptr, unsigned int) // 2{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, th; ROT2D(PARAM(2),PARAM(3),PARAM(4)) r=-( x2*x2 +x2*y2 -2*x2*PARAM_Y -PARAM_X*y2 +y2 ); return( min(PARAM(1), max(PARAM(0)*r,-PARAM(1))) );}DBL f_algbr_cyl4(FPUContext *ctx, DBL *ptr, unsigned int) // 3{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, th; ROT2D(PARAM(2),PARAM(3),PARAM(4)) r=-( x2*x2 +y2*y2 +2*x2*y2 +3*x2*PARAM_Y -y2*PARAM_Y ); return( min(PARAM(1), max(PARAM(0)*r,-PARAM(1))) );}DBL f_bicorn(FPUContext *ctx, DBL *ptr, unsigned int) // 4{ DBL r, r2,x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r = PARAM(1)*PARAM(1); r2=(x2 + z2 + 2*PARAM(1)*PARAM_Y - r); r = ( y2*(r - (x2 + z2)) - r2*r2); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_bifolia(FPUContext *ctx, DBL *ptr, unsigned int) // 5{ DBL r, x2=PARAM_X*PARAM_X, z2=PARAM_Z*PARAM_Z; r=x2+PARAM_Y*PARAM_Y+z2; r= -(r*r - PARAM(1)*(x2 + z2)*PARAM_Y); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_blob(FPUContext *ctx, DBL *ptr, unsigned int) // 6{ DBL r,r2,y2,z2,temp; /* blob */ y2=PARAM_Y*PARAM_Y; z2=PARAM_Z*PARAM_Z; r2= PARAM(0)*0.5; r=((PARAM_X+r2)*(PARAM_X+r2)+y2+z2)*PARAM(2)*PARAM(2); if (r>1.) r=1.; temp=( (PARAM_X - r2)*(PARAM_X - r2)+y2+z2)*PARAM(4)*PARAM(4); if (temp>1) temp=1.; return(-(PARAM(1) *(1.-r) *(1.-r) +PARAM(3) *(1.-temp)*(1.-temp) ));}DBL f_blob2(FPUContext *ctx, DBL *ptr, unsigned int) // 7{ DBL r,x2,y2,z2; x2=PARAM_X*PARAM_X; y2=PARAM_Y*PARAM_Y; z2=PARAM_Z*PARAM_Z; /* f= f1 + f2 */ r=exp(-(x2+y2+z2)*PARAM(1)) + exp(-((PARAM_X-PARAM(0))*(PARAM_X-PARAM(0))+y2+z2)*PARAM(1)); return(PARAM(3)-r*PARAM(2));}DBL f_boy_surface(FPUContext *ctx, DBL *ptr, unsigned int) // 8{ DBL r, r2,ph,x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r2=1-PARAM_Z; ph=x2+y2; r= PARAM(1)*(64*r2*r2*r2*z2*PARAM_Z- 48*r2*r2*z2*(3*x2+3*y2+2*z2)+ 12*r2*PARAM_Z*(27*ph*ph-24*z2*ph+ 36*sqrt(2.0)*PARAM_Y*PARAM_Z*(y2-3*x2)+4*z2*z2)+ (9*x2+9*y2-2*z2)*(-81*ph*ph-72*z2*ph+ 108*sqrt(2.0)*PARAM_X*PARAM_Z*(x2-3*y2)+4*z2*z2) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_comma(FPUContext *ctx, DBL *ptr, unsigned int) // 9{ DBL r,th,temp; if ((PARAM_X==0)&&(PARAM_Z==0)) PARAM_X=0.0001; th=atan2(PARAM_Z,PARAM_X); r=sqrt( (PARAM_X + PARAM(0)*0.25)*(PARAM_X + PARAM(0)*0.25) + PARAM_Z*PARAM_Z); temp= cos(th*0.5) * PARAM(0)*0.5 -sqrt((r- PARAM(0)*0.75)*(r - PARAM(0)*0.75)+PARAM_Y*PARAM_Y); temp= min(PARAM_Z, temp); r= PARAM(0)*0.5 - sqrt( (PARAM_X - PARAM(0)*0.5)*(PARAM_X - PARAM(0)*0.5) + PARAM_Z*PARAM_Z + PARAM_Y*PARAM_Y) ; return(-(DBL)max(temp,r));}DBL f_cross_ellipsoids(FPUContext *ctx, DBL *ptr, unsigned int) // 10{ DBL r,x2,y2,z2; x2=PARAM_X*PARAM_X; y2=PARAM_Y*PARAM_Y; z2=PARAM_Z*PARAM_Z; /* f= max (f1, f2, f3) */ r=max(exp(-(x2*PARAM(0)+y2*PARAM(0)+z2)*PARAM(1)), exp(-(x2*PARAM(0)+y2+z2*PARAM(0))*PARAM(1))); r=max(r, exp(-(x2+y2*PARAM(0)+z2*PARAM(0))*PARAM(1))); return(PARAM(3)-r*PARAM(2));}DBL f_crossed_trough(FPUContext *ctx, DBL *ptr, unsigned int) // 11{ DBL r; r=( PARAM_X*PARAM_X * PARAM_Z*PARAM_Z - PARAM_Y ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_cubic_saddle(FPUContext *ctx, DBL *ptr, unsigned int) // 12{ DBL r; r=-(PARAM_X*PARAM_X*PARAM_X - PARAM_Y*PARAM_Y*PARAM_Y - PARAM_Z); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_cushion(FPUContext *ctx, DBL *ptr, unsigned int) // 13{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r= (z2*x2 - z2*z2 - 2*PARAM_Z*x2 + 2*z2*PARAM_Z + x2 - z2 -(x2 - PARAM_Z)*(x2 - PARAM_Z) - y2*y2 - 2*x2*y2 - y2*z2 + 2*y2*PARAM_Z + y2); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_devils_curve(FPUContext *ctx, DBL *ptr, unsigned int) // 14{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-(x2*x2 + 2*x2*z2 - 0.36*x2 - y2*y2 + 0.25*y2 + z2*z2); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_devils_curve_2d(FPUContext *ctx, DBL *ptr, unsigned int) // 15{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, th; ROT2D(PARAM(3),PARAM(4),PARAM(5)) r=-( x2 * (x2 - PARAM(1)*PARAM(1)) - y2 * (y2 - PARAM(2)*PARAM(2)) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_dupin_cyclid(FPUContext *ctx, DBL *ptr, unsigned int) // 16{ DBL r, r2,x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z, ph,th, p1,p2,p3,p4; r2=PARAM(5)*PARAM(5); ph=PARAM(2)*PARAM(2); th=PARAM(4)*PARAM(4); p1= ph - th - (PARAM(3) + PARAM(1))*(PARAM(3) + PARAM(1)); p2= ph - th - (PARAM(3) - PARAM(1))*(PARAM(3) - PARAM(1)); p3= PARAM(3)*PARAM(3); p4=PARAM(1)*PARAM(1); r=-( p1*p2* (x2*x2+y2*y2+z2*z2)+ 2*(p1*p2* (x2*y2+x2*z2+y2*z2))+ 2*r2*((-th-p3 + ph+p4)* (2*PARAM_X*PARAM(3)+2*PARAM_Y*PARAM(4)-r2)-4*PARAM(4)*ph*PARAM_Y)* (x2+y2+z2)+ 4*r2*r2*(PARAM(3)*PARAM_X+PARAM(4)*PARAM_Y) *(-r2+PARAM(4)*PARAM_Y+PARAM(3)*PARAM_X)+ 4*r2*r2*p4*y2+ r2*r2*r2*r2); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_ellipsoid(FPUContext *ctx, DBL *ptr, unsigned int) // 17{ /* sphere */ return(sqrt(PARAM_X*PARAM_X * PARAM(0)*PARAM(0) + PARAM_Y*PARAM_Y * PARAM(1)*PARAM(1) + PARAM_Z*PARAM_Z * PARAM(2)*PARAM(2))); }DBL f_enneper(FPUContext *ctx, DBL *ptr, unsigned int) // 18{ DBL r, r2,x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; if (fabs(PARAM_Z)<0.2) PARAM_Z=0.2; r =((y2-x2)/(2*PARAM_Z)+2*z2/9+2/3); r2=((y2-x2)/(4*PARAM_Z)-(1/4)*(x2+y2+(8/9)*z2)+2/9); r=-( r*r*r -6*r2*r2) ; return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_flange_cover(FPUContext *ctx, DBL *ptr, unsigned int) // 19{ DBL r,x2,y2,z2; x2=PARAM_X*PARAM_X; y2=PARAM_Y*PARAM_Y; z2=PARAM_Z*PARAM_Z; /* f= f1 + f2 + f3 */ r=exp(-(x2*PARAM(0)+y2*PARAM(0)+z2)*PARAM(1)) + exp(-(x2*PARAM(0)+y2+z2*PARAM(0))*PARAM(1)) + exp(-(x2+y2*PARAM(0)+z2*PARAM(0))*PARAM(1)); return(PARAM(3)-r*PARAM(2));}DBL f_folium_surface(FPUContext *ctx, DBL *ptr, unsigned int) // 20{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=((y2 + z2) * (1+(PARAM(2) - 4*PARAM(1))*PARAM_X)+x2*(1 + PARAM(2))); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_folium_surface_2d(FPUContext *ctx, DBL *ptr, unsigned int) // 21{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z, th; ROT2D(PARAM(3),PARAM(4),PARAM(5)) r = -( y2 * (1 + (PARAM(2) - 4*PARAM(1))*PARAM_X) + x2*(1 + PARAM(2))); return( min(10., max(PARAM(0)*r, -10.)) );}DBL f_glob(FPUContext *ctx, DBL *ptr, unsigned int) // 22{ DBL r, x2=PARAM_X*PARAM_X; r= ( 0.5*x2*x2*PARAM_X + 0.5*x2*x2 - (PARAM_Y*PARAM_Y + PARAM_Z*PARAM_Z) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_heart(FPUContext *ctx, DBL *ptr, unsigned int) // 23{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-( (2*x2+y2+z2-1)*(2*x2+y2+z2-1)*(2*x2+y2+z2-1)- 0.1*x2*z2*PARAM_Z-y2*z2*PARAM_Z ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_helical_torus(FPUContext *ctx, DBL *ptr, unsigned int) // 24{ DBL r,r2, temp,th,ph; r2=sqrt(PARAM_X*PARAM_X+PARAM_Z*PARAM_Z); if ((PARAM_X==0)&&(PARAM_Z==0)) PARAM_X=0.000001; ph=atan2(PARAM_X,PARAM_Z); temp=atan2(r2 - PARAM(0), PARAM_Y); r=sqrt((r2-PARAM(0))*(r2-PARAM(0))+PARAM_Y*PARAM_Y); th=fmod(ph*PARAM(1)+temp*PARAM(2), TWO_M_PI); if (th<0) th+=TWO_M_PI; temp=atan2((th-M_PI)*PARAM(8), r-PARAM(5)); temp=cos(temp*PARAM(7)+ph*PARAM(9))*PARAM(3)+PARAM(4); r2=PARAM(4)+PARAM(5)-PARAM(0)-r2; th=th - M_PI; temp=temp-min(sqrt((r-PARAM(5))*(r-PARAM(5))+th*th*PARAM(6)), (r+PARAM(5))); return(-max(r2,temp));}DBL f_helix1(FPUContext *ctx, DBL *ptr, unsigned int) // 25{ DBL r, r2, temp, th, ph, x2; /* helix * * p[0] : number of helix * * p[1] : frequency * * p[2] : minor radius * * p[3] : major radius * * p[4] : shape parameter * * p[5] : cross section * * p[5] = 1: circle * * = 2: diamond * * < 1: rectangle(rounded) * * p[6] : rotation angle for p[5]<1 */ r = sqrt(PARAM_X * PARAM_X + PARAM_Z * PARAM_Z); if ((PARAM_X == 0) && (PARAM_Z == 0)) PARAM_X = 0.000001; th = atan2(PARAM_Z, PARAM_X); th = fmod(th * PARAM(0) + PARAM_Y * PARAM(1) * PARAM(0), TWO_M_PI); if (th < 0) th += TWO_M_PI; PARAM_Z = (th - M_PI) / PARAM(4) / (PARAM(1) * PARAM(0)); PARAM_X = r - PARAM(3); if (PARAM(5) == 1) r2 = sqrt(PARAM_X * PARAM_X + PARAM_Z * PARAM_Z); else { if (PARAM(6) != 0) { th = cos(PARAM(6) * M_PI_180); ph = sin(PARAM(6) * M_PI_180); x2 = PARAM_X * th - PARAM_Z * ph; PARAM_Z = PARAM_X * ph + PARAM_Z * th; PARAM_X = x2; } if (PARAM(5) != 0) { temp = 2. / PARAM(5); r2 = pow((pow(fabs(PARAM_X), temp) + pow(fabs(PARAM_Z), temp)), PARAM(5) *.5); } else r2 = max(fabs(PARAM_X), fabs(PARAM_Z)); } return (-PARAM(2) + min((PARAM(3) + r), r2));}DBL f_helix2(FPUContext *ctx, DBL *ptr, unsigned int) // 26{ DBL th, ph, x2, z2, r2, temp; /* helical shape for (minor radius>major radius * * cross section p[5] same as NFUNCTION = 6 */ th = PARAM_Y * PARAM(1); ph = cos(th); th = sin(th); x2 = PARAM_X - PARAM(3) * ph; z2 = PARAM_Z - PARAM(3) * th; PARAM_X = x2 * ph + z2 * th; PARAM_Z = (-x2 * th + z2 * ph); if (PARAM(5) == 1) return (sqrt(PARAM_X * PARAM_X + PARAM_Z * PARAM_Z) - PARAM(2)); if (PARAM(5) != 0) { temp = 2. / PARAM(5); r2 = pow((pow(fabs(PARAM_X), temp) + pow(fabs(PARAM_Z), temp)), PARAM(5) *.5); } else r2 = max(fabs(PARAM_X), fabs(PARAM_Z)); return (r2 - PARAM(2));}DBL f_hex_x(FPUContext *ctx, DBL *ptr, unsigned int) // 27{ DBL x1,y1,x2,y2, th; x1=fabs(fmod(fabs(PARAM_X), sqrt(3.0))-sqrt(3.0)/2); y1=fabs(fmod(fabs(PARAM_Y), 3)-1.5); x2=sqrt(3.0)/2-x1; y2=1.5-y1; if ((x1*x1+y1*y1)>(x2*x2+y2*y2)) { x1=x2; y1=y2; } if ((x1==0)&&(y1==0)) PARAM_X=0.000001; th=atan2(y1,x1); if (th<M_PI/6) return(x1); else { x1=cos(M_PI/3)*x1+sin(M_PI/3)*y1; return(x1); }}DBL f_hex_y(FPUContext *ctx, DBL *ptr, unsigned int) // 28{ DBL x1,y1,x2,y2, th; x1=fabs(fmod(fabs(PARAM_X), sqrt(3.0))-sqrt(3.0)/2); y1=fabs(fmod(fabs(PARAM_Y), 3)-1.5); x2=sqrt(3.0)/2-x1; y2=1.5-y1; if ((x1*x1+y1*y1)>(x2*x2+y2*y2)) { x1=x2; y1=y2; } if ((x1==0)&&(y1==0)) PARAM_X=0.000001; th=atan2(y1,x1); if (th<M_PI/6) return(y1); else { y1=-sin(M_PI/3)*x1+cos(M_PI/3)*y1; return(fabs(y1)); }}DBL f_hetero_mf(FPUContext *ctx, DBL *ptr, unsigned int fn) // 29{ DBL signal; VECTOR V1; DBL rem; Make_Vector(V1, PARAM_X, PARAM_Y, PARAM_Z); int ngen = (int)PARAM(5) & 3; signal = (Noise(V1, ngen)*2.0 - 1.0) + PARAM(3); VScaleEq(V1, PARAM(1)); DBL p1_2_mp0 = pow(PARAM(1), -PARAM(0)), ea = p1_2_mp0; for (int i = 1; i < PARAM(2); i++) { // make a noisy increment and scale it by f^(-H) DBL inc = ((Noise(V1,ngen)*2.0 - 1.0) + PARAM(3)) * ea; // scale the increment by (the current signal)^PARAM(4) at V1 // so that PARAM(4)=0 gives 'straight' 1/f, PARAM(4)=1 // gives heterogenous fractal, etc for (int p = (int) PARAM(4); p > 0; --p) inc *= signal; signal += inc; // go to next 'octave' VScaleEq(V1, PARAM(1)); ea *= p1_2_mp0; } rem = PARAM(2) - (int) PARAM(2); if(rem != 0.0) { // do something with fraction part of octave DBL inc = ((Noise(V1,ngen)*2.0-1.0) + PARAM(3)) * ea; signal += rem * inc * signal; } return signal;}DBL f_hunt_surface(FPUContext *ctx, DBL *ptr, unsigned int) // 30{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-( 4*(x2+y2+z2-13)*(x2+y2+z2-13)*(x2+y2+z2-13) + 27*(3*x2+y2-4*z2-12)*(3*x2+y2-4*z2-12) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_hyperbolic_torus(FPUContext *ctx, DBL *ptr, unsigned int) // 31{ DBL r, ph,th,x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; ph=PARAM(1)*PARAM(1); th=PARAM(2)*PARAM(2); r=-( x2*x2 + 2*x2*y2 - 2*x2*z2 - 2*(ph+th)*x2 + y2*y2 - 2*y2*z2 + 2*(ph-th)*y2 + z2*z2 + 2*(ph+th)*z2 + (ph-th)*(ph-th)); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_isect_ellipsoids(FPUContext *ctx, DBL *ptr, unsigned int) // 32{ DBL r,x2,y2,z2; x2=PARAM_X*PARAM_X; y2=PARAM_Y*PARAM_Y; z2=PARAM_Z*PARAM_Z; /* f= min (f1, f2, f3) */ r=min(exp(-(x2*PARAM(0)+y2*PARAM(0)+z2)*PARAM(1)), exp(-(x2*PARAM(0)+y2+z2*PARAM(0))*PARAM(1))); r=min(r, exp(-(x2+y2*PARAM(0)+z2*PARAM(0))*PARAM(1))); return(PARAM(3)-r*PARAM(2));}DBL f_kampyle_of_eudoxus(FPUContext *ctx, DBL *ptr, unsigned int) // 33{ DBL r, x2=PARAM_X*PARAM_X; r=-( (PARAM_Y*PARAM_Y + PARAM_Z*PARAM_Z) - PARAM(2)*PARAM(2) * x2*x2 + PARAM(2)*PARAM(2) * PARAM(1)*PARAM(1) * x2 ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_kampyle_of_eudoxus_2d(FPUContext *ctx, DBL *ptr, unsigned int) // 34{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, th; ROT2D(PARAM(3),PARAM(4),PARAM(5)) r=-( y2 - PARAM(2)*PARAM(2) * x2*x2 + PARAM(2)*PARAM(2) * PARAM(1)*PARAM(2) * x2 ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_klein_bottle(FPUContext *ctx, DBL *ptr, unsigned int) // 35{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-( (x2+y2+z2+2*PARAM_Y-1)*((x2+y2+z2-2*PARAM_Y-1)*(x2+y2+z2-2*PARAM_Y-1)-8*z2)+ 16*PARAM_X*PARAM_Z*(x2+y2+z2-2*PARAM_Y-1) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_kummer_surface_v1(FPUContext *ctx, DBL *ptr, unsigned int) // 36{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-( x2*x2+y2*y2+z2*z2-x2-y2-z2-x2*y2-x2*z2-y2*z2+1 ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_kummer_surface_v2(FPUContext *ctx, DBL *ptr, unsigned int) // 37{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-( x2*x2+y2*y2+z2*z2+PARAM(1)*(x2+y2+z2)+PARAM(2)*(x2*y2+x2*z2+y2*z2)+PARAM(3)*PARAM_X*PARAM_Y*PARAM_Z-1 ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_lemniscate_of_gerono(FPUContext *ctx, DBL *ptr, unsigned int) // 38{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-( x2*x2 - x2 + y2 + z2 ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_lemniscate_of_gerono_2d(FPUContext *ctx, DBL *ptr, unsigned int) // 39{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, th; ROT2D(PARAM(3),PARAM(4),PARAM(5)) r=-( y2 - PARAM(2)*PARAM(2) * PARAM(1)*PARAM(1) * x2 + PARAM(2)*PARAM(2)*x2*x2 ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_mesh1(FPUContext *ctx, DBL *ptr, unsigned int) // 40{ DBL th, ph, r, r2, temp; th = M_PI / PARAM(0); ph = M_PI / PARAM(1); r = fmod(PARAM_X, PARAM(0) * 2); if (r < 0) r += PARAM(0) * 2; r = fabs(r - PARAM(0)) * PARAM(2); r2 = (PARAM_Y - cos(PARAM_Z * ph) * PARAM(3)) * PARAM(4); temp = -sqrt(r2 * r2 + r * r); r = fmod(PARAM_X - PARAM(0), PARAM(0) * 2); if (r < 0) r += PARAM(0) * 2; r = fabs(r - PARAM(0)) * PARAM(2); r2 = (PARAM_Y + cos(PARAM_Z * ph) * PARAM(3)) * PARAM(4); temp = max(-sqrt(r2 * r2 + r * r), temp); r = fmod(PARAM_Z, PARAM(1) * 2); if (r < 0) r += PARAM(1) * 2; r = fabs(r - PARAM(1)) * PARAM(2); r2 = (PARAM_Y + cos(PARAM_X * th) * PARAM(3)) * PARAM(4); temp = max(-sqrt(r2 * r2 + r * r), temp); r = fmod(PARAM_Z - PARAM(1), PARAM(1) * 2); if (r < 0) r += PARAM(1) * 2; r = fabs(r - PARAM(1)) * PARAM(2); r2 = (PARAM_Y - cos(PARAM_X * th) * PARAM(3)) * PARAM(4); return (-max(-sqrt(r2 * r2 + r * r), temp));}DBL f_mitre(FPUContext *ctx, DBL *ptr, unsigned int) // 41{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-( 4*x2*(x2 + y2 + z2) - y2*(1 - y2 - z2) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_nodal_cubic(FPUContext *ctx, DBL *ptr, unsigned int) // 42{ DBL r; r=-( PARAM_Y*PARAM_Y*PARAM_Y + PARAM_Z*PARAM_Z*PARAM_Z - 6*PARAM_Y*PARAM_Z ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_odd(FPUContext *ctx, DBL *ptr, unsigned int) // 43{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=( z2*x2 - z2*z2 - 2*PARAM_Z*x2 + 2*z2*PARAM_Z + x2 - z2 - (x2 - PARAM_Z)*(x2 - PARAM_Z) - y2*y2 - 2*y2*x2 - y2*z2 + 2*y2*PARAM_Z + y2 ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_ovals_of_cassini(FPUContext *ctx, DBL *ptr, unsigned int) // 44{ DBL r, r2,x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r2 = (x2 + y2 + z2 + PARAM(1)*PARAM(1)); r = -(r2*r2 - PARAM(3)*PARAM(1)*PARAM(1)*(x2 + z2) - PARAM(2)*PARAM(2) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_paraboloid(FPUContext *ctx, DBL *ptr, unsigned int) // 45{ DBL r; r=-( PARAM_X*PARAM_X - PARAM_Y + PARAM_Z*PARAM_Z); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_parabolic_torus(FPUContext *ctx, DBL *ptr, unsigned int) // 46{ DBL r, ph,th, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; ph=PARAM(1)*PARAM(1); th=PARAM(2)*PARAM(2); r=-( x2*x2 + 2*x2*y2 - 2*x2*PARAM_Z - (ph+th)*x2 + y2*y2 - 2*y2*PARAM_Z + (ph-th)*y2 + z2 + (ph+th)*PARAM_Z + (ph-th)* (ph-th) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_ph(FPUContext *ctx, DBL *ptr, unsigned int) // 47{ return( atan2(sqrt(PARAM_X*PARAM_X + PARAM_Z*PARAM_Z ),PARAM_Y) ); }DBL f_pillow(FPUContext *ctx, DBL *ptr, unsigned int) // 48{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=( x2*x2 + y2*y2 + z2*z2 - (x2 + y2 + z2) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_piriform(FPUContext *ctx, DBL *ptr, unsigned int) // 49{ DBL r, x2=PARAM_X*PARAM_X; r=-( (x2*x2 - x2*PARAM_X) + PARAM_Y*PARAM_Y + PARAM_Z*PARAM_Z ); return( -min(10., max(PARAM(0)*r,-10.)) );}DBL f_piriform_2d(FPUContext *ctx, DBL *ptr, unsigned int) // 50{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y,th; ROT2D(PARAM(4),PARAM(5),PARAM(6)) r=-( PARAM_Y*PARAM_Y - PARAM(1) * PARAM(3)* PARAM(3) * x2*PARAM_X - PARAM(2) * PARAM(3)*PARAM(3) * x2*x2 ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_poly4(FPUContext *ctx, DBL *ptr, unsigned int) // 51{ /* cylindrical shape */ DBL y2, temp; y2=PARAM_Y*PARAM_Y; temp=(PARAM(0)+PARAM(1)*PARAM_Y+PARAM(2)*y2 + PARAM(3)*y2*PARAM_Y +PARAM(4)*y2*y2); if (temp<-5.) temp=-5.; return (-temp+sqrt(PARAM_X*PARAM_X+PARAM_Z*PARAM_Z));}DBL f_polytubes(FPUContext *ctx, DBL *ptr, unsigned int) // 52{ DBL x2,y2,z2,r,r2; int i; y2=PARAM_Y*PARAM_Y; r2=PARAM(1)+PARAM(2)*PARAM_Y+PARAM(3)*y2+PARAM(4)*y2*PARAM_Y+ PARAM(5)*y2*y2; r=-10000; for (i=0;i<(int)PARAM(0);i++) { x2 = PARAM_X - r2*sin(2*M_PI/PARAM(0)*i); z2 = PARAM_Z - r2*cos(2*M_PI/PARAM(0)*i); r= max(r,-sqrt( x2*x2 + z2*z2 ) ); } return(-r);}DBL f_quantum(FPUContext *ctx, DBL *ptr, unsigned int) // 53{ DBL r, th,temp; /* well known function in quantum mechanics */ if ((PARAM_X==0)&&(PARAM_Z==0)) PARAM_X=1e-6; r=sqrt(PARAM_X*PARAM_X+PARAM_Z*PARAM_Z); th=cos(atan2(r,PARAM_Y)); r=sqrt(r*r+PARAM_Y*PARAM_Y)*2.; temp=r*r*exp(-r*0.33333333)*(3.*th*th-1); return (temp*temp-12.0)*(-0.1);}DBL f_quartic_paraboloid(FPUContext *ctx, DBL *ptr, unsigned int) // 54{ DBL r, x2=PARAM_X*PARAM_X, z2=PARAM_Z*PARAM_Z; r=-( x2*x2 + z2*z2 - PARAM_Y ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_quartic_saddle(FPUContext *ctx, DBL *ptr, unsigned int) // 55{ DBL r, x2=PARAM_X*PARAM_X, z2=PARAM_Z*PARAM_Z; r=-( x2*x2 - z2*z2 - PARAM_Y ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_quartic_cylinder(FPUContext *ctx, DBL *ptr, unsigned int) // 56{ DBL r, x2=PARAM_X*PARAM_X, z2=PARAM_Z*PARAM_Z; r=-( (x2 + z2) * PARAM_Y*PARAM_Y + PARAM(2)*PARAM(2) * (x2 + z2) - PARAM(2)*PARAM(2) * PARAM(1)*PARAM(1) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_r(FPUContext *ctx, DBL *ptr, unsigned int) // 57{ return( sqrt(PARAM_X*PARAM_X + PARAM_Y*PARAM_Y + PARAM_Z*PARAM_Z ) ); }DBL f_ridge(FPUContext *ctx, DBL *ptr, unsigned int) // 58{ VECTOR EPoint = { PARAM_X, PARAM_Y, PARAM_Z }; int i; DBL Lambda, Omega, l, o; DBL value, v, resid, tot = 1.0, off, ridge, scale, rscale; VECTOR temp; int Octaves; int ngen = (int)PARAM(5) & 0x03; Lambda = l = PARAM(0); Octaves = (int)PARAM(1); Omega = o = PARAM(2); off = PARAM(3); ridge = PARAM(4); rscale = 1.0 / max(ridge, 1.0 - ridge); scale = 1.0 / max(off, 1.0 - off); resid = PARAM(1) - (DBL)Octaves; v = fabs(Noise(EPoint,ngen) - ridge) * rscale; value = (v - off); for (i = 2; i <= Octaves; i++) { VScale(temp,EPoint,l); v = fabs(Noise(temp,ngen) - ridge) * rscale; value += o * (v - off); tot += o; l *= Lambda; o *= Omega; } if (0.0 != resid) { VScale(temp,EPoint,l); v = fabs(Noise(temp,ngen) - ridge) * rscale; value += o * (v - off) * resid; tot += o * resid; } return (value * scale / tot);}DBL f_ridged_mf(FPUContext *ctx, DBL *ptr, unsigned int fn) // 59{ FunctionCode *f = ctx->functionvm->GetFunction(fn); DBL *ea,freq,signal,weight,result; int i; VECTOR V1; int ngen = (int)PARAM(5) & 0x03; V1[X]=PARAM_X; V1[Y]=PARAM_Y; V1[Z]=PARAM_Z; if (f->private_data == NULL) { ea = (DBL*)POV_MALLOC((PARAM(2) + 1)*sizeof(DBL), "exponent array"); freq = 1.0; for (i=0; i<=PARAM(2);i++) { ea[i]= pow(freq,-PARAM(0)); freq *= PARAM(1); } f->private_data = (void *)ea; } else { ea = (DBL *)(f->private_data); } signal = Noise(V1,ngen)*2.0-1.0; if (signal < 0.0 ) signal = -signal; signal = PARAM(3) - signal; signal *= signal; result = signal; weight = 1.0; for (i=1; i<PARAM(2); i++) { V1[X] *= PARAM(1); V1[Y] *= PARAM(1); V1[Z] *= PARAM(1); weight = signal * PARAM(4); if (weight > 1.0) weight = 1.0; if (weight < 0.0) weight = 0.0; signal = Noise(V1,ngen)*2.0-1.0; if (signal < 0.0 ) signal = -signal; signal = PARAM(3) - signal; signal *= signal; signal *= weight; result += signal * ea[i]; } return (result);}DBL f_rounded_box(FPUContext *ctx, DBL *ptr, unsigned int) // 60{ DBL x2, y2, z2, x3, y3, z3; x2 = PARAM(1) - PARAM(0); y2 = PARAM(2) - PARAM(0); z2 = PARAM(3) - PARAM(0); x3 = (PARAM_X < x2) ? 0 : (PARAM_X - x2); y3 = (PARAM_Y < y2) ? 0 : (PARAM_Y - y2); z3 = (PARAM_Z < z2) ? 0 : (PARAM_Z - z2); x2 = PARAM(0) - PARAM(1); y2 = PARAM(0) - PARAM(2); z2 = PARAM(0) - PARAM(3); PARAM_X = (PARAM_X > x2) ? 0 : (x2 - PARAM_X); PARAM_Y = (PARAM_Y > y2) ? 0 : (y2 - PARAM_Y); PARAM_Z = (PARAM_Z > z2) ? 0 : (z2 - PARAM_Z); PARAM_X = max(PARAM_X, x3); PARAM_Y = max(PARAM_Y, y3); PARAM_Z = max(PARAM_Z, z3); return (-PARAM(0) + sqrt(PARAM_X * PARAM_X + PARAM_Y * PARAM_Y + PARAM_Z * PARAM_Z) - 1e-6);}DBL f_sphere(FPUContext *ctx, DBL *ptr, unsigned int) // 61{ return (-PARAM(0) + sqrt(PARAM_X * PARAM_X + PARAM_Y * PARAM_Y + PARAM_Z * PARAM_Z));}DBL f_spikes(FPUContext *ctx, DBL *ptr, unsigned int) // 62{ DBL r,x2,y2,z2; x2=PARAM_X*PARAM_X; y2=PARAM_Y*PARAM_Y; z2=PARAM_Z*PARAM_Z; r= exp(-(x2+y2+z2)*PARAM(3))*PARAM(4) -exp(-(x2*PARAM(0)+y2*PARAM(0)+z2)*PARAM(1)) -exp(-(x2*PARAM(0)+y2+z2*PARAM(0))*PARAM(1)) -exp(-(x2+y2*PARAM(0)+z2*PARAM(0))*PARAM(1)); return(-r*PARAM(2));}DBL f_spikes_2d(FPUContext *ctx, DBL *ptr, unsigned int) // 63{ /* 2-D distribution */ return(-( PARAM(0)* cos(PARAM(1)*PARAM_X) * cos(PARAM(2)*PARAM_Z) * exp(- PARAM(3)*(PARAM_X*PARAM_X+PARAM_Z*PARAM_Z)) -PARAM_Y));}DBL f_spiral(FPUContext *ctx, DBL *ptr, unsigned int) // 64{ DBL r, r2, th, temp; /* spiral shape * * */ r = sqrt(PARAM_X * PARAM_X + PARAM_Z * PARAM_Z); if ((PARAM_X == 0) && (PARAM_Z == 0)) PARAM_X = 0.000001; th = atan2(PARAM_Z, PARAM_X); r = r + PARAM(0) * th / TWO_M_PI; r2 = fmod(r, PARAM(0)) - PARAM(0) * 0.5; if (PARAM(5) == 1) r2 = sqrt(r2 * r2 + PARAM_Y * PARAM_Y); else if (PARAM(5) != 0) { temp = 2 / PARAM(5); r2 = pow((pow(fabs(r2), temp) + pow(fabs(PARAM_Y), temp)), 1. / temp); } else r2 = max(fabs(r2), fabs(PARAM_Y)); r = sqrt(PARAM_X * PARAM_X + PARAM_Y * PARAM_Y + PARAM_Z * PARAM_Z); return (-min(PARAM(2) - r, PARAM(1) - min(r2, r)));}DBL f_steiners_roman(FPUContext *ctx, DBL *ptr, unsigned int) // 65{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-( x2*y2 + x2*z2 + y2*z2 + PARAM_X*PARAM_Y*PARAM_Z ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_strophoid(FPUContext *ctx, DBL *ptr, unsigned int) // 66{ DBL r, r2,x2=PARAM_X*PARAM_X; r2=PARAM(3)*PARAM(3); r=-((PARAM(2) - PARAM_X)*(PARAM_Y*PARAM_Y + PARAM_Z*PARAM_Z) - r2*PARAM(1)*x2 - r2*x2*PARAM_X ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_strophoid_2d(FPUContext *ctx, DBL *ptr, unsigned int) // 67{ DBL r, r2,x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, th; ROT2D(PARAM(4),PARAM(5),PARAM(6)) r2=PARAM(3)*PARAM(3); r=-( (PARAM(2) - PARAM_X)*y2 - r2*PARAM(1)*x2 - r2*x2*PARAM_X ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_superellipsoid(FPUContext *ctx, DBL *ptr, unsigned int) // 68{ DBL p = 2 / PARAM(0), n = 1 / PARAM(1); return 1 - pow((pow((pow(fabs(PARAM_X), p) + pow(fabs(PARAM_Y), p)), PARAM(0) * n) + pow(fabs(PARAM_Z), 2 * n)), PARAM(1) *.5);}DBL f_th(FPUContext *ctx, DBL *ptr, unsigned int) // 69{ return( atan2(PARAM_X,PARAM_Z) ); }DBL f_torus(FPUContext *ctx, DBL *ptr, unsigned int) // 70{ PARAM_X = sqrt(PARAM_X * PARAM_X + PARAM_Z * PARAM_Z) - PARAM(0); return -PARAM(1) + sqrt(PARAM_X * PARAM_X + PARAM_Y * PARAM_Y);}DBL f_torus2(FPUContext *ctx, DBL *ptr, unsigned int) // 71{ DBL r, ph,th,x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; ph=PARAM(1)*PARAM(1); th=PARAM(2)*PARAM(2); r=-( x2*x2 + y2*y2 + z2*z2 + 2*x2*y2 + 2*x2*z2 + 2*y2*z2 -2* (ph + th)* x2 + 2* (ph - th)* y2 -2* (ph + th)* z2 + (ph - th)*(ph - th) ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_torus_gumdrop(FPUContext *ctx, DBL *ptr, unsigned int) // 72{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, z2=PARAM_Z*PARAM_Z; r=-( 4*(x2*x2 + (y2 + z2)*(y2 + z2)) + 17 * x2 * (y2 + z2) - 20 * (x2 + y2 + z2) + 17); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_umbrella(FPUContext *ctx, DBL *ptr, unsigned int) // 73{ DBL r; r=-( PARAM_X*PARAM_X - PARAM_Y*PARAM_Z*PARAM_Z ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_witch_of_agnesi(FPUContext *ctx, DBL *ptr, unsigned int) // 74{ DBL r; r=-( PARAM(1) * (PARAM_Y - 1) + (PARAM_X*PARAM_X + PARAM_Z*PARAM_Z) * PARAM_Y ); return( min(10., max(PARAM(0)*r,-10.)) );}DBL f_witch_of_agnesi_2d(FPUContext *ctx, DBL *ptr, unsigned int) // 75{ DBL r, x2=PARAM_X*PARAM_X, y2=PARAM_Y*PARAM_Y, th; ROT2D(PARAM(3),PARAM(4),PARAM(5)) r=-( PARAM(1)*PARAM(1) * PARAM_Y + x2 * PARAM_Y - PARAM(2) ); return( min(10.0, max(PARAM(0)*r,-10.0)) );}DBL f_noise3d(FPUContext *ctx, DBL *ptr, unsigned int) // 76{ VECTOR Vec = { PARAM_X, PARAM_Y, PARAM_Z }; return Noise(Vec, ctx->threaddata->GetSceneData()->noiseGenerator);}DBL f_pattern(FPUContext *ctx, DBL *ptr, unsigned int fn) // 77{ VECTOR TPoint; VECTOR Vec = { PARAM_X, PARAM_Y, PARAM_Z }; FunctionCode *f = ctx->functionvm->GetFunction(fn); if(f->private_data == NULL) return 0.0; Warp_EPoint (TPoint, Vec, (TPATTERN *)(f->private_data)); return Evaluate_TPat((TPATTERN *)(f->private_data), TPoint, NULL, ctx->threaddata);}DBL f_noise_generator(FPUContext *ctx, DBL *ptr, unsigned int) // 78{ VECTOR Vec = { PARAM_X, PARAM_Y, PARAM_Z }; int ngen = (int)PARAM(0) & 0x03; return Noise(Vec, ngen);}DBL f_mandelbulb(FPUContext *ctx, DBL *ptr, unsigned int) // 79{ // PARAM_X .. PARAM_Z - x,y,z function parameters // PARAM(1) - power of n // PARAM(2) - iterations DBL x,y,z,newx, newy, newz, r, theta, phi, rn, thetan, phin, sinthetan, cx,cy,cz; int n,cnt,iterations; n = (int)PARAM(0); iterations = PARAM(1); cnt = iterations; cx = PARAM_X; cy = PARAM_Y; cz = PARAM_Z; x = 0.0; y = 0.0; z = 0.0; while(cnt > 0 && x*x + y*y + z*z < 2.0) { r = sqrt(x*x + y*y + z*z); theta = atan2(sqrt(x*x + y*y), z); phi = atan2(y,x); rn = pow(r,n); thetan = theta * n; phin = phi * n; sinthetan = sin(thetan); newx = cx + rn * sinthetan * (cos(phin)); newy = cy + rn * sinthetan * (sin(phin)); newz = cz + rn * (cos(thetan)); // now we reassign them x = newx; y = newy; z = newz; cnt--; } if (cnt <= 0) { return -1.0; } else {return 0.0; } /* return (1.0*cnt) / iterations; */}void f_pigment(FPUContext *ctx, DBL *ptr, unsigned int fn, unsigned int sp) // 0{ VECTOR Vec = { PARAM_N_X(5), PARAM_N_Y(5), PARAM_N_Z(5) }; COLOUR Col; FunctionCode *f = ctx->functionvm->GetFunction(fn); if(f->private_data == NULL) { ctx->SetLocal(sp + pRED, 0.0); ctx->SetLocal(sp + pGREEN, 0.0); ctx->SetLocal(sp + pBLUE, 0.0); ctx->SetLocal(sp + pFILTER, 0.0); ctx->SetLocal(sp + pTRANSM, 0.0); return; } Compute_Pigment(Col, (PIGMENT *)(f->private_data), Vec, NULL, ctx->threaddata); ctx->SetLocal(sp + pRED, Col[pRED]); ctx->SetLocal(sp + pGREEN, Col[pGREEN]); ctx->SetLocal(sp + pBLUE, Col[pBLUE]); ctx->SetLocal(sp + pFILTER, Col[pFILTER]); ctx->SetLocal(sp + pTRANSM, Col[pTRANSM]);}void f_transform(FPUContext *ctx, DBL *ptr, unsigned int fn, unsigned int sp) // 1{ VECTOR Vec = { PARAM_N_X(3), PARAM_N_Y(3), PARAM_N_Z(3) }; VECTOR Result; FunctionCode *f = ctx->functionvm->GetFunction(fn); if(f->private_data == NULL) { ctx->SetLocal(sp + X, 0.0); ctx->SetLocal(sp + Y, 0.0); ctx->SetLocal(sp + Z, 0.0); return; } MTransPoint(Result, Vec, (TRANSFORM *)(f->private_data)); ctx->SetLocal(sp + X, Result[X]); ctx->SetLocal(sp + Y, Result[Y]); ctx->SetLocal(sp + Z, Result[Z]);}void f_spline(FPUContext *ctx, DBL *ptr, unsigned int fn, unsigned int sp) // 2{ EXPRESS Result; FunctionCode *f = ctx->functionvm->GetFunction(fn); int Terms; if(f->private_data == NULL) { ctx->SetLocal(sp + X, 0.0); ctx->SetLocal(sp + Y, 0.0); ctx->SetLocal(sp + Z, 0.0); return; } Terms = ((SPLINE *)(f->private_data))->Terms; Get_Spline_Val((SPLINE *)(f->private_data), PARAM_N_X(Terms), Result, &Terms); ctx->SetLocal(sp + X, Result[X]); ctx->SetLocal(sp + Y, Result[Y]); if(Terms > 2) { ctx->SetLocal(sp + Z, Result[Z]); if(Terms > 3) { ctx->SetLocal(sp + T, Result[T]); if(Terms > 4) ctx->SetLocal(sp + T + 1, Result[T + 1]); } }}}