diff --git a/build/make.linux64.test.py b/build/make.linux64.test.py
index 5328ebe..0c2e6d1 100644
--- a/build/make.linux64.test.py
+++ b/build/make.linux64.test.py
@@ -17,7 +17,7 @@ builddir = "./build_linux64"
 doinstall = True   #copies the build_output to the install dir when finished
 cc = "nvcc" #compiler
 cflags = "-dc --compiler-options '-fPIC -O3'"
-libraries = "-l{} -lamsculib3.linux64 -lpthread".format(libname)
+libraries = "-l{} -lamsculib3.linux64 -lamsstring4.linux64 -lamsmathutil25.linux64 -lpthread".format(libname)
 libdirs = "-L{} -L{}/lib -L{}/lib".format(builddir,commondir,depdir)
 linkerflags = " -Xlinker=-rpath,."
 srcexts = [".c",".cpp",".cu"]
diff --git a/build_linux64/default.inp b/build_linux64/default.inp
new file mode 100644
index 0000000..9072cea
--- /dev/null
+++ b/build_linux64/default.inp
@@ -0,0 +1,33 @@
+#A default input file
+
+//comments are either c or python style
+
+avariable = bla
+
+anothervar = 1.2,1.3, 1.5 //this is a comment
+
+operatingmode = 1
+//1 - STL output
+//2 - count triangles
+//3 - count area
+
+verbosity = 0
+
+//domain
+xyz_min = -1.2,-1.2,-1.2
+xyz_max = 1.2,1.2,1.2
+Nx = 32
+Ny = 32
+Nz = 32
+
+//execution parallelism pars
+nblockdiv = 32
+nthreads = 256
+
+//function specific parameters
+//not all pars apply to all functions
+
+isjulia = 0
+juliac = 0.0, 0.0, 0.0
+exponent = 8.0
+
diff --git a/build_linux64/libamscudafractallevelset.linux64.a b/build_linux64/libamscudafractallevelset.linux64.a
index 81c9471..930c1c9 100644
Binary files a/build_linux64/libamscudafractallevelset.linux64.a and b/build_linux64/libamscudafractallevelset.linux64.a differ
diff --git a/build_linux64/objstore/amscufracls_bagoftriangles.o b/build_linux64/objstore/amscufracls_bagoftriangles.o
index 2ea4133..f9ab147 100644
Binary files a/build_linux64/objstore/amscufracls_bagoftriangles.o and b/build_linux64/objstore/amscufracls_bagoftriangles.o differ
diff --git a/build_linux64/objstore/amscufracls_function.o b/build_linux64/objstore/amscufracls_function.o
index f2b2311..e7f3de4 100644
Binary files a/build_linux64/objstore/amscufracls_function.o and b/build_linux64/objstore/amscufracls_function.o differ
diff --git a/build_linux64/objstore/amscufracls_marchingtets.o b/build_linux64/objstore/amscufracls_marchingtets.o
index fe7a4d2..d205be5 100644
Binary files a/build_linux64/objstore/amscufracls_marchingtets.o and b/build_linux64/objstore/amscufracls_marchingtets.o differ
diff --git a/build_linux64/objstore/amscufracls_misc.o b/build_linux64/objstore/amscufracls_misc.o
index 98586e8..9b34f28 100644
Binary files a/build_linux64/objstore/amscufracls_misc.o and b/build_linux64/objstore/amscufracls_misc.o differ
diff --git a/build_linux64/objstore/amscufracls_mtcube.o b/build_linux64/objstore/amscufracls_mtcube.o
index 5600fc2..d2807bc 100644
Binary files a/build_linux64/objstore/amscufracls_mtcube.o and b/build_linux64/objstore/amscufracls_mtcube.o differ
diff --git a/build_linux64/objstore/amscufracls_parseinputfile.o b/build_linux64/objstore/amscufracls_parseinputfile.o
new file mode 100644
index 0000000..a77cf24
Binary files /dev/null and b/build_linux64/objstore/amscufracls_parseinputfile.o differ
diff --git a/build_linux64/objstore/amscufracls_template.o b/build_linux64/objstore/amscufracls_template.o
index d93ce9f..c0862fc 100644
Binary files a/build_linux64/objstore/amscufracls_template.o and b/build_linux64/objstore/amscufracls_template.o differ
diff --git a/build_linux64/objstore/amscufracls_tetrahedron.o b/build_linux64/objstore/amscufracls_tetrahedron.o
index fe640ac..4ce75a7 100644
Binary files a/build_linux64/objstore/amscufracls_tetrahedron.o and b/build_linux64/objstore/amscufracls_tetrahedron.o differ
diff --git a/build_linux64/test b/build_linux64/test
index e341a0b..97a5741 100644
Binary files a/build_linux64/test and b/build_linux64/test differ
diff --git a/include/amscudafractallevelset/amscfls_functiondefs.cuh b/include/amscudafractallevelset/amscfls_functiondefs.cuh
new file mode 100644
index 0000000..f6deff0
--- /dev/null
+++ b/include/amscudafractallevelset/amscfls_functiondefs.cuh
@@ -0,0 +1,14 @@
+#ifndef __AMSCLFLS_FUNCTIONDEFS_CUH__
+#define __AMSCLFLS_FUNCTIONDEFS_CUH__
+
+#ifdef AMSCLS_MANDELBROT
+struct main_functionpars
+{
+    public:
+    int fractaliter;
+    float threshold;
+    float freezethresh;
+};
+#endif
+
+#endif
\ No newline at end of file
diff --git a/include/amscudafractallevelset/amscfls_parseinput.cuh b/include/amscudafractallevelset/amscfls_parseinput.cuh
new file mode 100644
index 0000000..6150770
--- /dev/null
+++ b/include/amscudafractallevelset/amscfls_parseinput.cuh
@@ -0,0 +1,43 @@
+#ifndef __AMSCFLS_PARSEINPUT_CUH__
+#define __AMSCFLS_PARSEINPUT_CUH__
+
+namespace amscuda
+{
+namespace fractallevelset
+{
+
+    struct inputpars
+    {
+        public:
+        cuvec3f xyz_min;
+        cuvec3f xyz_max;
+        int Nx;
+        int Ny;
+        int Nz;
+        int nblockdiv;
+        int nthreads;
+
+        bool isjulia;
+        cuvec3f juliac;
+        float exponent;
+
+        
+
+        int operatingmode;
+        int verbosity;
+        
+        __host__ __device__ inputpars();
+    };
+
+    int parse_input(const char *fname, inputpars *pars);
+
+    void inputmain(int argc, char *argv[]);
+
+    
+
+
+
+};//end namespaces
+};
+
+#endif
\ No newline at end of file
diff --git a/include/amscudafractallevelset/amscudafractallevelset.cuh b/include/amscudafractallevelset/amscudafractallevelset.cuh
index 9a112c8..f8a344d 100644
--- a/include/amscudafractallevelset/amscudafractallevelset.cuh
+++ b/include/amscudafractallevelset/amscudafractallevelset.cuh
@@ -38,6 +38,11 @@ class cuvec4f;
 // #define AMSCU_CONST
 // #endif
 
+#define AMSCLS_MANDELBROT
+//#define AMSCLS_JULIA
+//#define AMSCLS_PPKFRACTAL
+
+
 namespace amscuda
 {
 namespace fractallevelset
@@ -45,18 +50,10 @@ namespace fractallevelset
 
 __device__ __host__ int isignf(float f);
 
-struct main_functionpars
-{
-    public:
-    int fractaliter;
-    float threshold;
-    float freezethresh;
-};
+#include <amscudafractallevelset/amscfls_functiondefs.cuh>
 
 __host__ __device__ float main_function(cuvec3f pos, main_functionpars &pars);
 
-
-
 class bagoftriangles
 {
     public:
@@ -160,18 +157,26 @@ struct marchingtet_pars
     __device__ __host__ marchingtet_pars();
 };
 
+}; //end namespaces
+};
 
+#include <amscudafractallevelset/amscfls_parseinput.cuh>
+
+namespace amscuda
+{
+namespace fractallevelset
+{
 ///////////
 // Tests //
 ///////////
 void test_tetorientation();
 void test_mtcube_triangulation();
 void test_marchingtets();
+void hausdorff_test();
 
 
 
-
-};
+};//end namespaces
 };
 
 #endif
diff --git a/make_linux.py b/make_linux.py
index 84a1f5c..03ef5a6 100644
--- a/make_linux.py
+++ b/make_linux.py
@@ -18,5 +18,5 @@ for o in obj_list:
     os.remove('{}'.format(o))
 
 os.chdir('./build_linux64')
-callproc('./test')
+callproc('./test ./default.inp')
 os.chdir('..')
diff --git a/src/amscudafractallevelset/amscufracls_function.cu b/src/amscudafractallevelset/amscufracls_function.cu
index 2099c1f..9688b8d 100644
--- a/src/amscudafractallevelset/amscufracls_function.cu
+++ b/src/amscudafractallevelset/amscufracls_function.cu
@@ -6,16 +6,61 @@ namespace fractallevelset
 {
 
 
+// __host__ __device__ float main_function(cuvec3f pos, main_functionpars &pars)
+// {
+//     float ret = 0.0f;
+
+//     ret = pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]-pars.threshold;
+
+//     return -ret;
+// };
+
+#ifdef AMSCLS_MANDELBROT
+__host__ __device__ float main_function(cuvec3f pos, main_functionpars &pars)
+{
+    float ret = 0.0f;
+    float wr,wtheta,wphi;
+    cuvec3f w;
+    int I;
+
+    w = pos;
+
+    for(I=0;I<pars.fractaliter;I++)
+    {
+        wr = cuvec3f_norm(w);
+        wtheta = ::acos(w.y/wr);
+        wphi = ::atan2(w.x,w.z);
+        if(wr>pars.freezethresh) break;
+        wr = ::pow(wr,8.0f);
+        wtheta = wtheta*8.0f;
+        wphi = wphi*8.0f;
+        w = cuvec3f(wr*sin(wtheta)*sin(wphi),wr*cos(wtheta),wr*sin(wtheta)*cos(wphi));
+        w = w + pos;
+    }
+
+    ret = cuvec3f_norm(w)-pars.threshold;
+
+    return ret;
+};
+#endif
+
+#ifdef AMSCLS_JULIA
 __host__ __device__ float main_function(cuvec3f pos, main_functionpars &pars)
 {
     float ret = 0.0f;
 
-    ret = pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]-pars.threshold;
-
-    return -ret;
+    return ret;
 };
+#endif
 
+#ifdef AMSCLS_PPKFRACTAL
+__host__ __device__ float main_function(cuvec3f pos, main_functionpars &pars)
+{
+    float ret = 0.0f;
 
+    return ret;
+};
+#endif
     
 };//end namespace fractallevelset
 };// end namespace amscuda
diff --git a/src/amscudafractallevelset/amscufracls_marchingtets.cu b/src/amscudafractallevelset/amscufracls_marchingtets.cu
index 54e17ba..282584a 100644
--- a/src/amscudafractallevelset/amscufracls_marchingtets.cu
+++ b/src/amscudafractallevelset/amscufracls_marchingtets.cu
@@ -298,7 +298,7 @@ namespace fractallevelset
         }
         
         marchingtets_gpu_areacount<<<nblocks,nthreads>>>(mtpars,dareabuffer);
-        printf("debug: launching marchingtets_gpu_areacount with <%d,%d>\n",nblocks,nthreads);
+        //printf("debug: launching marchingtets_gpu_areacount with <%d,%d>\n",nblocks,nthreads);
         cudaError_t err = cudaSuccess;
         cudaDeviceSynchronize();
         err = cudaGetLastError();
@@ -355,11 +355,11 @@ namespace fractallevelset
 
     void test_marchingtets()
     {
-        const char *fnameout = "./testout.stl";
+        const char *fnameout = "./testout2.stl";
         bagoftriangles mesh;
         main_functionpars fpars;
         fpars.fractaliter = 64;
-        fpars.threshold = 0.5f;
+        fpars.threshold = 2.0f;
         fpars.freezethresh = 10.0f;
 
         marchingtet_pars mtpars;
@@ -389,5 +389,64 @@ namespace fractallevelset
 
     }
     
+
+
+    void hausdorff_test1a(int fractaliter, int N, float &deltat, float &area, float &narea)
+    {
+        main_functionpars fpars;
+        fpars.fractaliter = fractaliter;
+        fpars.threshold = 2.0f;
+        fpars.freezethresh = 10.0f;
+
+        marchingtet_pars mtpars;
+        mtpars.xyz_min = cuvec3f(-1.2,-1.2,-1.2);
+        mtpars.xyz_max = cuvec3f(1.2,1.2,1.2);
+        mtpars.Nx = N;
+        mtpars.Ny = N;
+        mtpars.Nz = N;
+    
+        mtpars.pars = fpars;
+
+        float t0,t1;
+
+        float dx = (mtpars.xyz_max[0]-mtpars.xyz_min[0])/((float)N);
+
+        //printf("debug: dA=%1.3f\n",dx*dx);
+
+        t0 = util::time_msec();
+        area = marchingtets_area_gpu(mtpars);
+        t1 = util::time_msec();
+        narea = area/dx/dx;
+
+        deltat=t1-t0;
+
+        // printf("marchingtets_area_gpu: %1.3f msec runtime %d psize\n",t1-t0,mtpars.Nx*mtpars.Ny*mtpars.Nz);
+        // printf("area=%1.3f, area/4pi=%1.3f\n",area,area/4.0f/pif);
+    }
+
+    void hausdorff_test()
+    {
+        int fiters[] = {32,64,128,256};
+        //int fiter;
+        int N = 32;
+
+        int I,J;
+
+        float deltat,area,narea;
+
+        for(I=0;I<4;I++)
+        {
+            N=32;
+            for(J=0;J<5;J++)
+            {
+                hausdorff_test1a(fiters[I],N,deltat,area,narea);
+                printf("%d,%d,%1.6g,%1.6g,%1.6g\n",fiters[I],N,deltat,area,narea);
+                N*=2;
+            }
+            printf("\n");
+        }
+
+    }
+
 };//end namespace fractallevelset
 };// end namespace amscuda
diff --git a/src/amscudafractallevelset/amscufracls_parseinputfile.cu b/src/amscudafractallevelset/amscufracls_parseinputfile.cu
new file mode 100644
index 0000000..00bc33c
--- /dev/null
+++ b/src/amscudafractallevelset/amscufracls_parseinputfile.cu
@@ -0,0 +1,236 @@
+#include <amscudafractallevelset/amscudafractallevelset.cuh>
+#include <amsmathutil25/amsmathutil25.hpp>
+#include <amsstring4/amsstring4.hpp>
+
+using namespace ams;
+
+namespace amscuda
+{
+namespace fractallevelset
+{
+
+    static const int omode_stloutput = 1;
+    static const int omode_counttriangles = 2;
+    static const int omode_countarea = 3;
+    
+    
+
+    inputpars::inputpars()
+    {
+        xyz_min = cuvec3f(-1.2,-1.2,-1.2);
+        xyz_max = cuvec3f(1.2,1.2,1.2);
+
+        Nx = 32;
+        Ny = 32;
+        Nz = 32;
+
+        nblockdiv = 32;
+        nthreads = 256;
+
+        isjulia = 0;
+        juliac = cuvec3f(0,0,0);
+        exponent = 8.0f;
+
+        operatingmode = omode_stloutput;
+        verbosity = 0; //0 silent, 1
+
+
+    }
+
+    
+
+    int load_inputfile(const char *fname, amsarray<amsarray<amsstring>> *stringpars)
+    {
+        int ret = 1;
+        int res;
+        amsstring ftxt,ln;
+        amsarray<amsstring> lns;
+        amsarray<amsstring> pair;
+        int I;
+
+        if(stringpars==NULL)
+        {
+            printf("parse_input error: stringpars is NULL\n");
+            return -1;
+        }
+
+        stringpars->resize(0);
+
+        res = readtextfile(fname,&ftxt);
+        if(res<0)
+        {
+            printf("parse_input error: could not open %s for reading.\n",fname);
+            return -1;
+        }
+
+        lns = ftxt.splitlines();
+        for(I=0;I<lns.length;I++)
+        {
+            ln = lns[I].splitfirst("#")[0];
+            ln = ln.splitfirst("//")[0];
+            if(ln.length>0)
+            {
+                pair = ln.splitfirst("=");
+                if(pair[0].length>0 && pair[1].length>0)
+                {
+                    pair[0] = pair[0].stripwhitespace();
+                    pair[1] = pair[1].stripwhitespace();
+                    stringpars->append(pair);
+                }
+            }
+        }
+
+        // printf("debug: parfile pars:\n");
+        // for(I=0;I<stringpars->length;I++)
+        // {
+        //     printf("%d:\t'%s'\t'%s'\n",I,stringpars->at(I)[0].cstring,stringpars->at(I)[1].cstring);
+        // }
+
+        return ret;
+    }
+
+    cuvec3f parse_as_cuvec3f(amsstring s)
+    {
+        int I;
+        cuvec3f ret;
+        amsarray<amsstring> alphanums;
+        
+        alphanums = s.splitalphanum();
+        for(I=0;I<alphanums.length&&I<3;I++)
+        {
+            ret[I] = (float)(alphanums[I].strtonum());
+        }
+
+        return ret;
+    }
+
+    void parsearch_cuvec3f(amsstring varname, amsarray<amsarray<amsstring>> *stringpars, cuvec3f *val)
+    {
+        int I;
+        for(I=0;I<stringpars->length;I++)
+        {
+            if(varname==stringpars->at(I)[0])
+            {
+                *val = parse_as_cuvec3f(stringpars->at(I)[1]);
+            }
+        }
+        return;
+    }
+
+    void parsearch_float(amsstring varname, amsarray<amsarray<amsstring>> *stringpars, float *val)
+    {
+        int I;
+        for(I=0;I<stringpars->length;I++)
+        {
+            if(varname==stringpars->at(I)[0])
+            {
+                *val = (float) stringpars->at(I)[1].strtonum();
+            }
+        }
+        return;
+    }
+
+    void parsearch_int(amsstring varname, amsarray<amsarray<amsstring>> *stringpars, int *val)
+    {
+        int I;
+        for(I=0;I<stringpars->length;I++)
+        {
+            if(varname==stringpars->at(I)[0])
+            {
+                *val = (int) stringpars->at(I)[1].strtonum();
+            }
+        }
+        return;
+    }
+
+    void parsearch_bool(amsstring varname, amsarray<amsarray<amsstring>> *stringpars, bool *val)
+    {
+        int I;
+        for(I=0;I<stringpars->length;I++)
+        {
+            if(varname==stringpars->at(I)[0])
+            {
+                *val = (bool) stringpars->at(I)[1].strtonum();
+            }
+        }
+        return;
+    }
+
+    int parse_input(const char *fname, inputpars *pars)
+    {
+        int ret = 1;
+        int res;
+        amsarray<amsarray<amsstring>> stringpars;
+        res = load_inputfile(fname,&stringpars);
+
+        parsearch_cuvec3f("xyz_min",&stringpars,&(pars->xyz_min));
+        parsearch_cuvec3f("xyz_max",&stringpars,&(pars->xyz_max));
+        parsearch_int("Nx",&stringpars,&(pars->Nx));
+        parsearch_int("Ny",&stringpars,&(pars->Ny));
+        parsearch_int("Nz",&stringpars,&(pars->Nz));
+        parsearch_int("nblockdiv",&stringpars,&(pars->nblockdiv));
+        parsearch_int("nthreads",&stringpars,&(pars->nthreads));
+        parsearch_bool("juliac",&stringpars,&(pars->isjulia));
+        parsearch_int("operatingmode",&stringpars,&(pars->operatingmode));
+        parsearch_int("verbosity",&stringpars,&(pars->verbosity));
+
+
+        return ret;
+    }
+
+    void inputmain(int argc, char *argv[])
+    {
+        int res;
+        amsstring fnameout;
+        inputpars pars;
+        marchingtet_pars mtpars;
+        main_functionpars mfpars;
+
+        if(argc<3)
+        {
+            fnameout = "./output.stl";
+        }
+        else
+        {
+            fnameout = argv[2];
+        }
+
+        if(argc<2)
+        {
+            
+        }
+        else
+        {
+            res = parse_input(argv[1],&pars);
+            if(res<0) return;
+        }
+
+        //Fill operatingstructs
+        mtpars.xyz_min = pars.xyz_min;
+        mtpars.xyz_max = pars.xyz_max;
+        mtpars.Nx = pars.Nx;
+        mtpars.Ny = pars.Ny;
+        mtpars.Nz = pars.Nz;
+        mtpars.pars = mfpars;
+
+        //execute model based on operatingmode
+        if(pars.operatingmode==omode_stloutput)
+        {
+
+        }
+        else if(pars.operatingmode==omode_counttriangles)
+        {
+
+        }
+        else if(pars.operatingmode==omode_countarea)
+        {
+
+        }
+
+
+    }
+
+
+    
+};//end namespace fractallevelset
+};// end namespace amscuda
diff --git a/src/main.cu b/src/main.cu
index 5f0e292..5d153ef 100644
--- a/src/main.cu
+++ b/src/main.cu
@@ -12,8 +12,10 @@ int main(int argc, char* argv[])
     //test_tetorientation();
     //test_mtcube_triangulation();
 
-    test_marchingtets();
+    //test_marchingtets();
+    //hausdorff_test();
 
+    inputmain(argc,argv);
 
     return 0;
 }
\ No newline at end of file