diff --git a/python/frep-C.py b/python/frep-C.py
new file mode 100644
index 0000000000000000000000000000000000000000..8999ac5e2b47d66944ee358a96828e3ae287f04e
--- /dev/null
+++ b/python/frep-C.py
@@ -0,0 +1,132 @@
+#!/usr/bin/env python3
+#
+# frep-C.py
+#    functional representation to C++ solver
+#
+# usage:
+#    pcb.py | frep-C.py [dpi [filename]]
+#
+# Neil Gershenfeld 8/23/21
+# (c) Massachusetts Institute of Technology 2021
+#
+# This work may be reproduced, modified, distributed,
+# performed, and displayed for any purpose, but must
+# acknowledge this project. Copyright is retained and
+# must be preserved. The work is provided as is; no
+# warranty is provided, and users accept all liability.
+#
+
+#
+# import
+#
+import json,sys,os
+from numpy import *
+from PIL import Image
+#
+# read input
+#
+frep = json.load(sys.stdin)
+#
+# check arguments
+#
+if (frep['type'] != 'RGB'):
+   print('types other than RGB not (yet) supported')
+   sys.exit()
+if (len(sys.argv) == 1):
+   print('output to out.png at 300 DPI')
+   filename = 'out.png'
+   dpi = 300
+elif (len(sys.argv) == 2):
+   dpi = sys.argv[1]
+   filename = 'out.png'
+   print('output to out.png at '+dpi+'DPI')
+   dpi = int(dpi)
+elif (len(sys.argv) == 3):
+   dpi = sys.argv[1]
+   filename = sys.argv[2]
+   print('output to '+filename+' at '+dpi+' DPI')
+   dpi = int(dpi)
+#
+# define variables
+#
+xmin = frep['xmin']
+xmax = frep['xmax']
+ymin = frep['ymin']
+ymax = frep['ymax']
+units = float(frep['mm_per_unit'])
+delta = (25.4/dpi)/units
+fn = frep['function']
+fn = fn.replace('math.','')
+#
+# write program
+#
+file = open("frep-C.cpp",'w')
+file.write(
+f"""
+#include <iostream>
+#include <cmath>
+#include <png.h>
+using namespace std;
+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 units = {units};
+   int dpi = {dpi};
+   float delta = (25.4/dpi)/units;
+   int nx = (xmax-xmin)/delta;
+   int ny = (ymax-ymin)/delta;
+   int *m = (int*) malloc(nx*ny*sizeof(int));
+   for (int ix = 0; ix < nx; ++ix) {{
+      float x = xmin+ix*delta;
+      for (int iy = 0; iy < ny; ++iy) {{
+         float y = ymin+iy*delta;
+         m[iy*nx+ix] = fn(x,y,0);
+         }}
+      }}
+   FILE *file;
+   file = fopen("{filename}","wb");
+	png_structp pngfile;
+	png_infop pnginfo;
+	png_bytep pngrow;
+	pngfile = png_create_write_struct(PNG_LIBPNG_VER_STRING,NULL,NULL,NULL);
+	pnginfo = png_create_info_struct(pngfile);
+	png_init_io(pngfile,file);
+	png_set_IHDR(pngfile,pnginfo,nx,ny,
+	   8,PNG_COLOR_TYPE_RGBA,PNG_INTERLACE_NONE,
+	   PNG_COMPRESSION_TYPE_BASE,PNG_FILTER_TYPE_BASE);
+   png_set_pHYs(pngfile,pnginfo,{1000*dpi/25.4},{1000*dpi/25.4},PNG_RESOLUTION_METER);
+	png_write_info(pngfile,pnginfo);
+	pngrow = (png_bytep) malloc(4*nx*sizeof(png_byte));
+	for (int iy = (ny-1); iy >= 0; --iy) {{
+		for (int ix = 0; ix < nx; ++ix) {{
+         pngrow[4*ix] = m[iy*nx+ix] & 255;
+         pngrow[4*ix+1] = (m[iy*nx+ix] >> 8) & 255; 
+         pngrow[4*ix+2] = (m[iy*nx+ix] >> 16) & 255;
+         pngrow[4*ix+3] = 255;
+         }}
+		png_write_row(pngfile,pngrow);
+	   }}
+	png_write_end(pngfile,NULL);
+	fclose(file);
+   }}
+""")
+file.close()
+#
+# compile
+#
+print("compile ...")
+os.system("time g++ frep-C.cpp -o frep-C -lm -lpng -O -ffast-math")
+#
+# execute
+#
+print("execute ...")
+os.system("time ./frep-C")
+#
+# clean up
+#
+os.system("rm -f frep-C.cpp frep-C")