diff --git a/python/frep-C.py b/python/frep-C.py
index ce21cc38aec5380f8fb3d09ea0a657dbfff8333f..d57f7b0c92e8e3d711b778e5ffa6ff65c9e1ee82 100755
--- a/python/frep-C.py
+++ b/python/frep-C.py
@@ -1,13 +1,12 @@
 #!/usr/bin/env python3
 #
 # frep-C.py
-#    functional representation to C++ solver
+#    functional representation to C++ multi-threaded solver
 #
 # usage:
 #    pcb.py | frep-C.py [dpi [filename]]
 #
-# Neil Gershenfeld 8/23/21
-# (c) Massachusetts Institute of Technology 2021
+# Neil Gershenfeld 11/26/22
 #
 # This work may be reproduced, modified, distributed,
 # performed, and displayed for any purpose, but must
@@ -71,34 +70,43 @@ file.write(
 f"""
 #include <iostream>
 #include <cmath>
+#include <thread>
 #include <png.h>
+//
 using namespace std;
+//
+float xmin = {xmin};
+float xmax = {xmax};
+float ymin = {ymin};
+float ymax = {ymax};
+float zmin = {zmin};
+float zmax = {zmax};
+float units = {units};
+int dpi = {dpi};
+float delta = (25.4/dpi)/units;
+int nx = (xmax-xmin)/delta;
+int ny = (ymax-ymin)/delta;
+int *m = (int*) calloc(nx*ny,sizeof(int));
+float layers[] = {{{layers}}};
+int nthreads = std::thread::hardware_concurrency();
+//
 int fn(float X,float Y,float Z) {{
    return ({fn});
    }}
-int main(int argc, char** argv) {{
-   float xmin = {xmin};
-   float xmax = {xmax};
-   float ymin = {ymin};
-   float ymax = {ymax};
-   float zmin = {zmin};
-   float zmax = {zmax};
-   float units = {units};
-   int dpi = {dpi};
-   float delta = (25.4/dpi)/units;
-   int nx = (xmax-xmin)/delta;
-   int ny = (ymax-ymin)/delta;
-   int *m = (int*) calloc(nx*ny,sizeof(int));
-   float layers[] = {{{layers}}};
+//
+void calc(int nx,int ny,int nthreads,int thread) {{
    int intensity;
-   for (int layer = 0; layer < {nlayers}; ++layer) {{
+   for (int layer = 0; layer < 2; ++layer) {{
       float z = layers[layer];
-      cout << "   z =" << z << endl;
+      if (thread == 0)
+         cout << "   z = " << z << endl;
       if (zmin == zmax)
-         intensity = 0xFFFFFF;
+         intensity = 0xffffff;
       else
          intensity = ((int) (255*(z-zmin)/(zmax-zmin))) | (255 << 8) | (255 << 16);
-      for (int iy = 0; iy < ny; ++iy) {{
+      int iystart = thread*ny/nthreads;
+      int iyend = (thread+1)*ny/nthreads;
+      for (int iy = iystart; iy < iyend; ++iy) {{
          float y = ymin+iy*delta;
          for (int ix = 0; ix < nx; ++ix) {{
             float x = xmin+ix*delta;
@@ -106,6 +114,16 @@ int main(int argc, char** argv) {{
             }}
          }}
       }}
+   }}
+//
+int main(int argc, char** argv) {{
+   cout << "   calculate " << nx << "x" << ny << " with " << nthreads << " threads" << endl;
+   std::thread threads[nthreads];
+   for (int i = 0; i < nthreads; ++i)
+      threads[i] = std::thread(calc,nx,ny,nthreads,i);
+   for (int i = 0; i < nthreads; ++i)
+      threads[i].join();
+   //
    FILE *file;
    file = fopen("{filename}","wb");
 	png_structp pngfile;
@@ -138,7 +156,7 @@ file.close()
 # compile
 #
 print("compile ...")
-os.system("time g++ frep-C.cpp -o frep-C -lm -lpng -O -ffast-math")
+os.system("time g++ frep-C.cpp -o frep-C -lm -lpng -O -ffast-math -pthread")
 #
 # execute
 #
@@ -147,4 +165,4 @@ os.system("time ./frep-C")
 #
 # clean up
 #
-os.system("rm -f frep-C.cpp frep-C")
+#os.system("rm -f frep-C.cpp frep-C")