diff --git a/.github/workflows/build-linux-clazy.yml b/.github/workflows/build-linux-clazy.yml
index cac8dd19..e4bed7ed 100644
--- a/.github/workflows/build-linux-clazy.yml
+++ b/.github/workflows/build-linux-clazy.yml
@@ -37,6 +37,6 @@ jobs:
       # all level 1 checks but ignore clazy-no-connect-by-name
       # no parallel make with clazy to not mess log
       - run: |
-          export CLAZY_IGNORE_DIRS=".*usr.*|.*bezier.*|.*boost.*|.*SingleApplication.*|.*spdlog.*|.*zernike.*" \
+          export CLAZY_IGNORE_DIRS=".*usr.*|.*bezier.*|.*boost.*|.*SingleApplication.*|.*spdlog.*|.*zernike.*|.*cnpy.*" \
           && export CLAZY_CHECKS="level1,no-connect-by-name,function-args-by-value,function-args-by-ref,incorrect-emit,old-style-connect" \
           && make
diff --git a/.github/workflows/build-linux.yml b/.github/workflows/build-linux.yml
index 75d55b65..b483b24f 100644
--- a/.github/workflows/build-linux.yml
+++ b/.github/workflows/build-linux.yml
@@ -49,4 +49,4 @@ jobs:
           tidy-checks: '' # rely solely on .clang-tidy file
           tidy-review: true
           passive-reviews: true
-          ignore: 'bezier|boost|SingleApplication|spdlog|zernike|moc_*|ui_*|qwt*'
+          ignore: 'bezier|boost|SingleApplication|spdlog|zernike|cnpy|moc_*|ui_*|qwt*'
diff --git a/DFTFringe.pro b/DFTFringe.pro
index 4adcc4e2..e94e9cab 100644
--- a/DFTFringe.pro
+++ b/DFTFringe.pro
@@ -56,6 +56,7 @@ win32 {
     LIBS += -L$$PWD\..\build_openCV\install\x64\mingw\bin -llibopencv_imgcodecs4120
     LIBS += -L$$PWD\..\build_openCV\install\x64\mingw\bin -llibopencv_imgproc4120
     LIBS += -ldbghelp # for SetUnhandledExceptionFilter
+    LIBS += -lz       # zip compression library needed for cnpy.cpp
 
 
     # This is for armadillo to not use wrapper. See https://gitlab.com/conradsnicta/armadillo-code#6-linux-and-macos-compiling-and-linking
@@ -81,6 +82,7 @@ unix: !mac {
     LIBS += -lopencv_imgcodecs
     LIBS += -lopencv_imgproc
     LIBS += -L/usr/local/qwt-6.3.0/lib -lqwt
+    LIBS += -lz       # zip compression library needed for cnpy.cpp
 }
 
 # MAC ##############
@@ -143,7 +145,7 @@ RESOURCES += DFTResources.qrc
 
 TRANSLATIONS += dftfringe_fr.ts
 
-INCLUDEPATH += ./bezier ./SingleApplication ./zernike
+INCLUDEPATH += ./bezier ./SingleApplication ./zernike ./cnpy
 
 SOURCES += SingleApplication/singleapplication.cpp \
     SingleApplication/singleapplication_p.cpp \
@@ -164,6 +166,7 @@ SOURCES += SingleApplication/singleapplication.cpp \
     ccswappeddlg.cpp \
     circlefit.cpp \
     circleoutline.cpp \
+    cnpy/cnpy.cpp \
     colorchannel.cpp \
     colorchanneldisplay.cpp \
     colormapviewerdlg.cpp \
@@ -285,6 +288,7 @@ HEADERS += bezier/bezier.h \
     circle.h \
     circleoutline.h \
     circleutils.h \
+    cnpy/cnpy.h \
     colorchannel.h \
     colorchanneldisplay.h \
     colormapviewerdlg.h \
diff --git a/DFTFringe_Dale.pro b/DFTFringe_Dale.pro
index 24b5d0e7..76fd813f 100644
--- a/DFTFringe_Dale.pro
+++ b/DFTFringe_Dale.pro
@@ -39,6 +39,7 @@ SOURCES += main.cpp \
     mainwindow.cpp \
     igramarea.cpp \
     circleoutline.cpp \
+    cnpy/cnpy.cpp \
     graphicsutilities.cpp \
     dfttools.cpp \
     dftarea.cpp \
@@ -156,6 +157,7 @@ HEADERS  += mainwindow.h \
     edgeplot.h \
     IgramArea.h \
     circleoutline.h \
+    cnpy/cnpy.h \
     graphicsutilities.h \
     dfttools.h \
     dftarea.h \
@@ -266,7 +268,7 @@ HEADERS  += mainwindow.h \
     SingleApplication/singleapplication.h \
     SingleApplication/singleapplication_p.h
 
-INCLUDEPATH += ./bezier ./SingleApplication ./zernike
+INCLUDEPATH += ./bezier ./SingleApplication ./zernike ./cnpy
 
 FORMS    += mainwindow.ui \
     annulushelpdlg.ui \
@@ -378,6 +380,7 @@ LIBS += D:\lapack\build64\bin\liblapack.dll
 LIBS += D:\lapack\build64\bin\libblas.dll
 
 LIBS += -ldbghelp # for SetUnhandledExceptionFilter
+LIBS += -lz       # zip compression library needed for cnpy.cpp
 }
 
 
diff --git a/DFTFringe_QT5.pro b/DFTFringe_QT5.pro
index 79b3336e..65573323 100644
--- a/DFTFringe_QT5.pro
+++ b/DFTFringe_QT5.pro
@@ -56,7 +56,7 @@ win32 {
     LIBS += -L..\build_openCV\install\x64\mingw\bin -llibopencv_imgcodecs460
     LIBS += -L..\build_openCV\install\x64\mingw\bin -llibopencv_imgproc460
     LIBS += -ldbghelp # for SetUnhandledExceptionFilter
-
+    LIBS += -lz       # zip compression library needed for cnpy.cpp
 
     # This is for armadillo to not use wrapper. See https://gitlab.com/conradsnicta/armadillo-code#6-linux-and-macos-compiling-and-linking
     DEFINES += ARMA_DONT_USE_WRAPPER
@@ -80,6 +80,8 @@ unix: !mac {
     LIBS += -lopencv_imgproc
     LIBS += -lopencv_imgproc
     LIBS += -lqwt-qt5
+    LIBS += -lz       # zip compression library needed for cnpy.cpp
+
 }
 
 # MAC ##############
@@ -142,7 +144,7 @@ RESOURCES += DFTResources.qrc
 
 TRANSLATIONS += dftfringe_fr.ts
 
-INCLUDEPATH += ./bezier ./SingleApplication ./zernike
+INCLUDEPATH += ./bezier ./SingleApplication ./zernike ./cnpy
 
 SOURCES += SingleApplication/singleapplication.cpp \
     SingleApplication/singleapplication_p.cpp \
@@ -163,6 +165,7 @@ SOURCES += SingleApplication/singleapplication.cpp \
     ccswappeddlg.cpp \
     circlefit.cpp \
     circleoutline.cpp \
+    cnpy/cnpy.cpp \
     colorchannel.cpp \
     colorchanneldisplay.cpp \
     colormapviewerdlg.cpp \
@@ -284,6 +287,7 @@ HEADERS += bezier/bezier.h \
     circle.h \
     circleoutline.h \
     circleutils.h \
+    cnpy/cnpy.h \
     colorchannel.h \
     colorchanneldisplay.h \
     colormapviewerdlg.h \
diff --git a/LICENSES/MIT.txt b/LICENSES/MIT.txt
index 0474e915..9003e035 100644
--- a/LICENSES/MIT.txt
+++ b/LICENSES/MIT.txt
@@ -15,6 +15,9 @@ Copyright (C) 2016 Gabi Melman.
 zernike (zapm.cpp)
 Copyright (C) 2022 Michael Peck
 
+cnpy (cnpy.h cnpy.cpp)
+Copyright (c) Carl Rogers, 2011
+
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
 in the Software without restriction, including without limitation the rights
diff --git a/RevisionHistory.html b/RevisionHistory.html
index a17bf02b..83944cbc 100644
--- a/RevisionHistory.html
+++ b/RevisionHistory.html
@@ -10,6 +10,14 @@
     <h1>DFTFringe Version History</h1>
     <ul>
 
+        <li>Version 8.4.0</li>
+        <ul>
+            <li>Ronchi and Foucault displays have new circular grid feature to measure zones better (right click on ronchi)</li>
+            <li>Ronchi compare now has blink mode that switches between 2 ronchis (first select at least 2 ronchis, then right click on the ronchi view, then select "compare" button, then "blink" button)</li>
+            <li>Can read Wavefront Pro format wavefront files (.npz files)</li>
+            <li>Fixed bug in subtract wavefront feature introduced in version 7.4.0 where sometimes it subtracts the wrong wavefront (but correctly showed you what it did on the right side).  Bug happened on A-B if B is farther down in the wavefront list than A.</li>
+        </ul>
+
         <li>Version 8.3.2</li>
         <ul>
             <li>Fixed bug introduced in 8.3.1 in Ronchi Foucault where it was always doing autocollimate mode</li>
diff --git a/cnpy/LICENSE.txt b/cnpy/LICENSE.txt
new file mode 100644
index 00000000..e60eadbc
--- /dev/null
+++ b/cnpy/LICENSE.txt
@@ -0,0 +1,21 @@
+The MIT License
+
+Copyright (c) Carl Rogers, 2011
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
diff --git a/cnpy/cnpy.cpp b/cnpy/cnpy.cpp
new file mode 100644
index 00000000..f42275dd
--- /dev/null
+++ b/cnpy/cnpy.cpp
@@ -0,0 +1,339 @@
+//Copyright (C) 2011  Carl Rogers
+//Released under MIT License
+//license available in LICENSE file, or at http://www.opensource.org/licenses/mit-license.php
+
+#include"cnpy.h"
+#include<complex>
+#include<cstdlib>
+#include<algorithm>
+#include<cstring>
+#include<iomanip>
+#include<stdint.h>
+#include<stdexcept>
+#include <regex>
+
+char cnpy::BigEndianTest() {
+    int x = 1;
+    return (((char *)&x)[0]) ? '<' : '>';
+}
+
+char cnpy::map_type(const std::type_info& t)
+{
+    if(t == typeid(float) ) return 'f';
+    if(t == typeid(double) ) return 'f';
+    if(t == typeid(long double) ) return 'f';
+
+    if(t == typeid(int) ) return 'i';
+    if(t == typeid(char) ) return 'i';
+    if(t == typeid(short) ) return 'i';
+    if(t == typeid(long) ) return 'i';
+    if(t == typeid(long long) ) return 'i';
+
+    if(t == typeid(unsigned char) ) return 'u';
+    if(t == typeid(unsigned short) ) return 'u';
+    if(t == typeid(unsigned long) ) return 'u';
+    if(t == typeid(unsigned long long) ) return 'u';
+    if(t == typeid(unsigned int) ) return 'u';
+
+    if(t == typeid(bool) ) return 'b';
+
+    if(t == typeid(std::complex<float>) ) return 'c';
+    if(t == typeid(std::complex<double>) ) return 'c';
+    if(t == typeid(std::complex<long double>) ) return 'c';
+
+    else return '?';
+}
+
+template<> std::vector<char>& cnpy::operator+=(std::vector<char>& lhs, const std::string rhs) {
+    lhs.insert(lhs.end(),rhs.begin(),rhs.end());
+    return lhs;
+}
+
+template<> std::vector<char>& cnpy::operator+=(std::vector<char>& lhs, const char* rhs) {
+    //write in little endian
+    size_t len = strlen(rhs);
+    lhs.reserve(len);
+    for(size_t byte = 0; byte < len; byte++) {
+        lhs.push_back(rhs[byte]);
+    }
+    return lhs;
+}
+
+void cnpy::parse_npy_header(unsigned char* buffer,size_t& word_size, std::vector<size_t>& shape, bool& fortran_order) {
+    //std::string magic_string(buffer,6);
+    uint8_t major_version = *reinterpret_cast<uint8_t*>(buffer+6);
+    uint8_t minor_version = *reinterpret_cast<uint8_t*>(buffer+7);
+    uint16_t header_len = *reinterpret_cast<uint16_t*>(buffer+8);
+    std::string header(reinterpret_cast<char*>(buffer+9),header_len);
+
+    size_t loc1, loc2;
+
+    //fortran order
+    loc1 = header.find("fortran_order")+16;
+    fortran_order = (header.substr(loc1,4) == "True" ? true : false);
+
+    //shape
+    loc1 = header.find("(");
+    loc2 = header.find(")");
+
+    std::regex num_regex("[0-9][0-9]*");
+    std::smatch sm;
+    shape.clear();
+
+    std::string str_shape = header.substr(loc1+1,loc2-loc1-1);
+    while(std::regex_search(str_shape, sm, num_regex)) {
+        shape.push_back(std::stoi(sm[0].str()));
+        str_shape = sm.suffix().str();
+    }
+
+    //endian, word size, data type
+    //byte order code | stands for not applicable. 
+    //not sure when this applies except for byte array
+    loc1 = header.find("descr")+9;
+    bool littleEndian = (header[loc1] == '<' || header[loc1] == '|' ? true : false);
+    assert(littleEndian);
+
+    //char type = header[loc1+1];
+    //assert(type == map_type(T));
+
+    std::string str_ws = header.substr(loc1+2);
+    loc2 = str_ws.find("'");
+    word_size = atoi(str_ws.substr(0,loc2).c_str());
+}
+
+void cnpy::parse_npy_header(FILE* fp, size_t& word_size, std::vector<size_t>& shape, bool& fortran_order) {  
+    char buffer[256];
+    size_t res = fread(buffer,sizeof(char),11,fp);       
+    if(res != 11)
+        throw std::runtime_error("parse_npy_header: failed fread");
+    std::string header = fgets(buffer,256,fp);
+    assert(header[header.size()-1] == '\n');
+
+    size_t loc1, loc2;
+
+    //fortran order
+    loc1 = header.find("fortran_order");
+    if (loc1 == std::string::npos)
+        throw std::runtime_error("parse_npy_header: failed to find header keyword: 'fortran_order'");
+    loc1 += 16;
+    fortran_order = (header.substr(loc1,4) == "True" ? true : false);
+
+    //shape
+    loc1 = header.find("(");
+    loc2 = header.find(")");
+    if (loc1 == std::string::npos || loc2 == std::string::npos)
+        throw std::runtime_error("parse_npy_header: failed to find header keyword: '(' or ')'");
+
+    std::regex num_regex("[0-9][0-9]*");
+    std::smatch sm;
+    shape.clear();
+
+    std::string str_shape = header.substr(loc1+1,loc2-loc1-1);
+    while(std::regex_search(str_shape, sm, num_regex)) {
+        shape.push_back(std::stoi(sm[0].str()));
+        str_shape = sm.suffix().str();
+    }
+
+    //endian, word size, data type
+    //byte order code | stands for not applicable. 
+    //not sure when this applies except for byte array
+    loc1 = header.find("descr");
+    if (loc1 == std::string::npos)
+        throw std::runtime_error("parse_npy_header: failed to find header keyword: 'descr'");
+    loc1 += 9;
+    bool littleEndian = (header[loc1] == '<' || header[loc1] == '|' ? true : false);
+    assert(littleEndian);
+
+    //char type = header[loc1+1];
+    //assert(type == map_type(T));
+
+    std::string str_ws = header.substr(loc1+2);
+    loc2 = str_ws.find("'");
+    word_size = atoi(str_ws.substr(0,loc2).c_str());
+}
+
+void cnpy::parse_zip_footer(FILE* fp, uint16_t& nrecs, size_t& global_header_size, size_t& global_header_offset)
+{
+    std::vector<char> footer(22);
+    fseek(fp,-22,SEEK_END);
+    size_t res = fread(&footer[0],sizeof(char),22,fp);
+    if(res != 22)
+        throw std::runtime_error("parse_zip_footer: failed fread");
+
+    uint16_t disk_no, disk_start, nrecs_on_disk, comment_len;
+    disk_no = *(uint16_t*) &footer[4];
+    disk_start = *(uint16_t*) &footer[6];
+    nrecs_on_disk = *(uint16_t*) &footer[8];
+    nrecs = *(uint16_t*) &footer[10];
+    global_header_size = *(uint32_t*) &footer[12];
+    global_header_offset = *(uint32_t*) &footer[16];
+    comment_len = *(uint16_t*) &footer[20];
+
+    assert(disk_no == 0);
+    assert(disk_start == 0);
+    assert(nrecs_on_disk == nrecs);
+    assert(comment_len == 0);
+}
+
+cnpy::NpyArray load_the_npy_file(FILE* fp) {
+    std::vector<size_t> shape;
+    size_t word_size;
+    bool fortran_order;
+    cnpy::parse_npy_header(fp,word_size,shape,fortran_order);
+
+    cnpy::NpyArray arr(shape, word_size, fortran_order);
+    size_t nread = fread(arr.data<char>(),1,arr.num_bytes(),fp);
+    if(nread != arr.num_bytes())
+        throw std::runtime_error("load_the_npy_file: failed fread");
+    return arr;
+}
+
+cnpy::NpyArray load_the_npz_array(FILE* fp, uint32_t compr_bytes, uint32_t uncompr_bytes) {
+
+    std::vector<unsigned char> buffer_compr(compr_bytes);
+    std::vector<unsigned char> buffer_uncompr(uncompr_bytes);
+    size_t nread = fread(&buffer_compr[0],1,compr_bytes,fp);
+    if(nread != compr_bytes)
+        throw std::runtime_error("load_the_npy_file: failed fread");
+
+    int err;
+    z_stream d_stream;
+
+    d_stream.zalloc = Z_NULL;
+    d_stream.zfree = Z_NULL;
+    d_stream.opaque = Z_NULL;
+    d_stream.avail_in = 0;
+    d_stream.next_in = Z_NULL;
+    err = inflateInit2(&d_stream, -MAX_WBITS);
+
+    d_stream.avail_in = compr_bytes;
+    d_stream.next_in = &buffer_compr[0];
+    d_stream.avail_out = uncompr_bytes;
+    d_stream.next_out = &buffer_uncompr[0];
+
+    err = inflate(&d_stream, Z_FINISH);
+    err = inflateEnd(&d_stream);
+
+    std::vector<size_t> shape;
+    size_t word_size;
+    bool fortran_order;
+    cnpy::parse_npy_header(&buffer_uncompr[0],word_size,shape,fortran_order);
+
+    cnpy::NpyArray array(shape, word_size, fortran_order);
+
+    size_t offset = uncompr_bytes - array.num_bytes();
+    memcpy(array.data<unsigned char>(),&buffer_uncompr[0]+offset,array.num_bytes());
+
+    return array;
+}
+
+cnpy::npz_t cnpy::npz_load(std::string fname) {
+    FILE* fp = fopen(fname.c_str(),"rb");
+
+    if(!fp) {
+        throw std::runtime_error("npz_load: Error! Unable to open file "+fname+"!");
+    }
+
+    cnpy::npz_t arrays;  
+
+    while(1) {
+        std::vector<char> local_header(30);
+        size_t headerres = fread(&local_header[0],sizeof(char),30,fp);
+        if(headerres != 30)
+            throw std::runtime_error("npz_load: failed fread");
+
+        //if we've reached the global header, stop reading
+        if(local_header[2] != 0x03 || local_header[3] != 0x04) break;
+
+        //read in the variable name
+        uint16_t name_len = *(uint16_t*) &local_header[26];
+        std::string varname(name_len,' ');
+        size_t vname_res = fread(&varname[0],sizeof(char),name_len,fp);
+        if(vname_res != name_len)
+            throw std::runtime_error("npz_load: failed fread");
+
+        //erase the lagging .npy        
+        varname.erase(varname.end()-4,varname.end());
+
+        //read in the extra field
+        uint16_t extra_field_len = *(uint16_t*) &local_header[28];
+        if(extra_field_len > 0) {
+            std::vector<char> buff(extra_field_len);
+            size_t efield_res = fread(&buff[0],sizeof(char),extra_field_len,fp);
+            if(efield_res != extra_field_len)
+                throw std::runtime_error("npz_load: failed fread");
+        }
+
+        uint16_t compr_method = *reinterpret_cast<uint16_t*>(&local_header[0]+8);
+        uint32_t compr_bytes = *reinterpret_cast<uint32_t*>(&local_header[0]+18);
+        uint32_t uncompr_bytes = *reinterpret_cast<uint32_t*>(&local_header[0]+22);
+
+        if(compr_method == 0) {arrays[varname] = load_the_npy_file(fp);}
+        else {arrays[varname] = load_the_npz_array(fp,compr_bytes,uncompr_bytes);}
+    }
+
+    fclose(fp);
+    return arrays;  
+}
+
+cnpy::NpyArray cnpy::npz_load(std::string fname, std::string varname) {
+    FILE* fp = fopen(fname.c_str(),"rb");
+
+    if(!fp) throw std::runtime_error("npz_load: Unable to open file "+fname);
+
+    while(1) {
+        std::vector<char> local_header(30);
+        size_t header_res = fread(&local_header[0],sizeof(char),30,fp);
+        if(header_res != 30)
+            throw std::runtime_error("npz_load: failed fread");
+
+        //if we've reached the global header, stop reading
+        if(local_header[2] != 0x03 || local_header[3] != 0x04) break;
+
+        //read in the variable name
+        uint16_t name_len = *(uint16_t*) &local_header[26];
+        std::string vname(name_len,' ');
+        size_t vname_res = fread(&vname[0],sizeof(char),name_len,fp);      
+        if(vname_res != name_len)
+            throw std::runtime_error("npz_load: failed fread");
+        vname.erase(vname.end()-4,vname.end()); //erase the lagging .npy
+
+        //read in the extra field
+        uint16_t extra_field_len = *(uint16_t*) &local_header[28];
+        fseek(fp,extra_field_len,SEEK_CUR); //skip past the extra field
+        
+        uint16_t compr_method = *reinterpret_cast<uint16_t*>(&local_header[0]+8);
+        uint32_t compr_bytes = *reinterpret_cast<uint32_t*>(&local_header[0]+18);
+        uint32_t uncompr_bytes = *reinterpret_cast<uint32_t*>(&local_header[0]+22);
+
+        if(vname == varname) {
+            NpyArray array  = (compr_method == 0) ? load_the_npy_file(fp) : load_the_npz_array(fp,compr_bytes,uncompr_bytes);
+            fclose(fp);
+            return array;
+        }
+        else {
+            //skip past the data
+            uint32_t size = *(uint32_t*) &local_header[22];
+            fseek(fp,size,SEEK_CUR);
+        }
+    }
+
+    fclose(fp);
+
+    //if we get here, we haven't found the variable in the file
+    throw std::runtime_error("npz_load: Variable name "+varname+" not found in "+fname);
+}
+
+cnpy::NpyArray cnpy::npy_load(std::string fname) {
+
+    FILE* fp = fopen(fname.c_str(), "rb");
+
+    if(!fp) throw std::runtime_error("npy_load: Unable to open file "+fname);
+
+    NpyArray arr = load_the_npy_file(fp);
+
+    fclose(fp);
+    return arr;
+}
+
+
diff --git a/cnpy/cnpy.h b/cnpy/cnpy.h
new file mode 100644
index 00000000..4ea87a9c
--- /dev/null
+++ b/cnpy/cnpy.h
@@ -0,0 +1,269 @@
+//Copyright (C) 2011  Carl Rogers
+//Released under MIT License
+//license available in LICENSE file, or at http://www.opensource.org/licenses/mit-license.php
+
+#ifndef LIBCNPY_H_
+#define LIBCNPY_H_
+
+#include<string>
+#include<stdexcept>
+#include<sstream>
+#include<vector>
+#include<cstdio>
+#include<typeinfo>
+#include<iostream>
+#include<cassert>
+#include<zlib.h>
+#include<map>
+#include<memory>
+#include<stdint.h>
+#include<numeric>
+
+namespace cnpy {
+
+    struct NpyArray {
+        NpyArray(const std::vector<size_t>& _shape, size_t _word_size, bool _fortran_order) :
+            shape(_shape), word_size(_word_size), fortran_order(_fortran_order)
+        {
+            num_vals = 1;
+            for(size_t i = 0;i < shape.size();i++) num_vals *= shape[i];
+            data_holder = std::shared_ptr<std::vector<char>>(
+                new std::vector<char>(num_vals * word_size));
+        }
+
+        NpyArray() : shape(0), word_size(0), fortran_order(0), num_vals(0) { }
+
+        template<typename T>
+        T* data() {
+            return reinterpret_cast<T*>(&(*data_holder)[0]);
+        }
+
+        template<typename T>
+        const T* data() const {
+            return reinterpret_cast<T*>(&(*data_holder)[0]);
+        }
+
+        template<typename T>
+        std::vector<T> as_vec() const {
+            const T* p = data<T>();
+            return std::vector<T>(p, p+num_vals);
+        }
+
+        size_t num_bytes() const {
+            return data_holder->size();
+        }
+
+        std::shared_ptr<std::vector<char>> data_holder;
+        std::vector<size_t> shape;
+        size_t word_size;
+        bool fortran_order;
+        size_t num_vals;
+    };
+   
+    using npz_t = std::map<std::string, NpyArray>; 
+
+    char BigEndianTest();
+    char map_type(const std::type_info& t);
+    template<typename T> std::vector<char> create_npy_header(const std::vector<size_t>& shape);
+    void parse_npy_header(FILE* fp,size_t& word_size, std::vector<size_t>& shape, bool& fortran_order);
+    void parse_npy_header(unsigned char* buffer,size_t& word_size, std::vector<size_t>& shape, bool& fortran_order);
+    void parse_zip_footer(FILE* fp, uint16_t& nrecs, size_t& global_header_size, size_t& global_header_offset);
+    npz_t npz_load(std::string fname);
+    NpyArray npz_load(std::string fname, std::string varname);
+    NpyArray npy_load(std::string fname);
+
+    template<typename T> std::vector<char>& operator+=(std::vector<char>& lhs, const T rhs) {
+        //write in little endian
+        for(size_t byte = 0; byte < sizeof(T); byte++) {
+            char val = *((char*)&rhs+byte); 
+            lhs.push_back(val);
+        }
+        return lhs;
+    }
+
+    template<> std::vector<char>& operator+=(std::vector<char>& lhs, const std::string rhs);
+    template<> std::vector<char>& operator+=(std::vector<char>& lhs, const char* rhs);
+
+
+    template<typename T> void npy_save(std::string fname, const T* data, const std::vector<size_t> shape, std::string mode = "w") {
+        FILE* fp = NULL;
+        std::vector<size_t> true_data_shape; //if appending, the shape of existing + new data
+
+        if(mode == "a") fp = fopen(fname.c_str(),"r+b");
+
+        if(fp) {
+            //file exists. we need to append to it. read the header, modify the array size
+            size_t word_size;
+            bool fortran_order;
+            parse_npy_header(fp,word_size,true_data_shape,fortran_order);
+            assert(!fortran_order);
+
+            if(word_size != sizeof(T)) {
+                std::cout<<"libnpy error: "<<fname<<" has word size "<<word_size<<" but npy_save appending data sized "<<sizeof(T)<<"\n";
+                assert( word_size == sizeof(T) );
+            }
+            if(true_data_shape.size() != shape.size()) {
+                std::cout<<"libnpy error: npy_save attempting to append misdimensioned data to "<<fname<<"\n";
+                assert(true_data_shape.size() != shape.size());
+            }
+
+            for(size_t i = 1; i < shape.size(); i++) {
+                if(shape[i] != true_data_shape[i]) {
+                    std::cout<<"libnpy error: npy_save attempting to append misshaped data to "<<fname<<"\n";
+                    assert(shape[i] == true_data_shape[i]);
+                }
+            }
+            true_data_shape[0] += shape[0];
+        }
+        else {
+            fp = fopen(fname.c_str(),"wb");
+            true_data_shape = shape;
+        }
+
+        std::vector<char> header = create_npy_header<T>(true_data_shape);
+        size_t nels = std::accumulate(shape.begin(),shape.end(),1,std::multiplies<size_t>());
+
+        fseek(fp,0,SEEK_SET);
+        fwrite(&header[0],sizeof(char),header.size(),fp);
+        fseek(fp,0,SEEK_END);
+        fwrite(data,sizeof(T),nels,fp);
+        fclose(fp);
+    }
+
+    template<typename T> void npz_save(std::string zipname, std::string fname, const T* data, const std::vector<size_t>& shape, std::string mode = "w")
+    {
+        //first, append a .npy to the fname
+        fname += ".npy";
+
+        //now, on with the show
+        FILE* fp = NULL;
+        uint16_t nrecs = 0;
+        size_t global_header_offset = 0;
+        std::vector<char> global_header;
+
+        if(mode == "a") fp = fopen(zipname.c_str(),"r+b");
+
+        if(fp) {
+            //zip file exists. we need to add a new npy file to it.
+            //first read the footer. this gives us the offset and size of the global header
+            //then read and store the global header.
+            //below, we will write the the new data at the start of the global header then append the global header and footer below it
+            size_t global_header_size;
+            parse_zip_footer(fp,nrecs,global_header_size,global_header_offset);
+            fseek(fp,global_header_offset,SEEK_SET);
+            global_header.resize(global_header_size);
+            size_t res = fread(&global_header[0],sizeof(char),global_header_size,fp);
+            if(res != global_header_size){
+                throw std::runtime_error("npz_save: header read error while adding to existing zip");
+            }
+            fseek(fp,global_header_offset,SEEK_SET);
+        }
+        else {
+            fp = fopen(zipname.c_str(),"wb");
+        }
+
+        std::vector<char> npy_header = create_npy_header<T>(shape);
+
+        size_t nels = std::accumulate(shape.begin(),shape.end(),1,std::multiplies<size_t>());
+        size_t nbytes = nels*sizeof(T) + npy_header.size();
+
+        //get the CRC of the data to be added
+        uint32_t crc = crc32(0L,(uint8_t*)&npy_header[0],npy_header.size());
+        crc = crc32(crc,(uint8_t*)data,nels*sizeof(T));
+
+        //build the local header
+        std::vector<char> local_header;
+        local_header += "PK"; //first part of sig
+        local_header += (uint16_t) 0x0403; //second part of sig
+        local_header += (uint16_t) 20; //min version to extract
+        local_header += (uint16_t) 0; //general purpose bit flag
+        local_header += (uint16_t) 0; //compression method
+        local_header += (uint16_t) 0; //file last mod time
+        local_header += (uint16_t) 0;     //file last mod date
+        local_header += (uint32_t) crc; //crc
+        local_header += (uint32_t) nbytes; //compressed size
+        local_header += (uint32_t) nbytes; //uncompressed size
+        local_header += (uint16_t) fname.size(); //fname length
+        local_header += (uint16_t) 0; //extra field length
+        local_header += fname;
+
+        //build global header
+        global_header += "PK"; //first part of sig
+        global_header += (uint16_t) 0x0201; //second part of sig
+        global_header += (uint16_t) 20; //version made by
+        global_header.insert(global_header.end(),local_header.begin()+4,local_header.begin()+30);
+        global_header += (uint16_t) 0; //file comment length
+        global_header += (uint16_t) 0; //disk number where file starts
+        global_header += (uint16_t) 0; //internal file attributes
+        global_header += (uint32_t) 0; //external file attributes
+        global_header += (uint32_t) global_header_offset; //relative offset of local file header, since it begins where the global header used to begin
+        global_header += fname;
+
+        //build footer
+        std::vector<char> footer;
+        footer += "PK"; //first part of sig
+        footer += (uint16_t) 0x0605; //second part of sig
+        footer += (uint16_t) 0; //number of this disk
+        footer += (uint16_t) 0; //disk where footer starts
+        footer += (uint16_t) (nrecs+1); //number of records on this disk
+        footer += (uint16_t) (nrecs+1); //total number of records
+        footer += (uint32_t) global_header.size(); //nbytes of global headers
+        footer += (uint32_t) (global_header_offset + nbytes + local_header.size()); //offset of start of global headers, since global header now starts after newly written array
+        footer += (uint16_t) 0; //zip file comment length
+
+        //write everything
+        fwrite(&local_header[0],sizeof(char),local_header.size(),fp);
+        fwrite(&npy_header[0],sizeof(char),npy_header.size(),fp);
+        fwrite(data,sizeof(T),nels,fp);
+        fwrite(&global_header[0],sizeof(char),global_header.size(),fp);
+        fwrite(&footer[0],sizeof(char),footer.size(),fp);
+        fclose(fp);
+    }
+
+    template<typename T> void npy_save(std::string fname, const std::vector<T> data, std::string mode = "w") {
+        std::vector<size_t> shape;
+        shape.push_back(data.size());
+        npy_save(fname, &data[0], shape, mode);
+    }
+
+    template<typename T> void npz_save(std::string zipname, std::string fname, const std::vector<T> data, std::string mode = "w") {
+        std::vector<size_t> shape;
+        shape.push_back(data.size());
+        npz_save(zipname, fname, &data[0], shape, mode);
+    }
+
+    template<typename T> std::vector<char> create_npy_header(const std::vector<size_t>& shape) {  
+
+        std::vector<char> dict;
+        dict += "{'descr': '";
+        dict += BigEndianTest();
+        dict += map_type(typeid(T));
+        dict += std::to_string(sizeof(T));
+        dict += "', 'fortran_order': False, 'shape': (";
+        dict += std::to_string(shape[0]);
+        for(size_t i = 1;i < shape.size();i++) {
+            dict += ", ";
+            dict += std::to_string(shape[i]);
+        }
+        if(shape.size() == 1) dict += ",";
+        dict += "), }";
+        //pad with spaces so that preamble+dict is modulo 16 bytes. preamble is 10 bytes. dict needs to end with \n
+        int remainder = 16 - (10 + dict.size()) % 16;
+        dict.insert(dict.end(),remainder,' ');
+        dict.back() = '\n';
+
+        std::vector<char> header;
+        header += (char) 0x93;
+        header += "NUMPY";
+        header += (char) 0x01; //major version of numpy format
+        header += (char) 0x00; //minor version of numpy format
+        header += (uint16_t) dict.size();
+        header.insert(header.end(),dict.begin(),dict.end());
+
+        return header;
+    }
+
+
+}
+
+#endif
\ No newline at end of file
diff --git a/mainwindow.cpp b/mainwindow.cpp
index f0301941..863bc7bc 100644
--- a/mainwindow.cpp
+++ b/mainwindow.cpp
@@ -333,7 +333,7 @@ void MainWindow::openWaveFrontonInit(QStringList args){
         if (pd.wasCanceled())
             break;
 
-        if (arg.endsWith(".wft", Qt::CaseInsensitive)){
+        if (arg.endsWith(".wft", Qt::CaseInsensitive) || arg.endsWith(".npz",Qt::CaseInsensitive)){
             pd.setLabelText(arg);
             try {
             m_surfaceManager->loadWavefront(arg);
@@ -663,7 +663,7 @@ QStringList MainWindow::SelectWaveFrontFiles(){
 
     QFileDialog dialog(this, "load wave front file", lastPath, tr("wft(*.wft)"));
     dialog.setFileMode(QFileDialog::ExistingFiles);
-    dialog.setNameFilter(tr("wft (*.wft)"));
+    dialog.setNameFilter(tr("wavefront files (*.wft *.npz)"));
 
     if (dialog.exec()) {
         QStringList fileNames = dialog.selectedFiles();
diff --git a/surfacemanager.cpp b/surfacemanager.cpp
index 5a94b983..21382f54 100644
--- a/surfacemanager.cpp
+++ b/surfacemanager.cpp
@@ -75,6 +75,7 @@
 #include "ui_oglrendered.h"
 #include "astigpolargraph.h"
 #include "utils.h"
+#include "cnpy.h"
 
 cv::Mat theMask;
 cv::Mat deb;
@@ -1133,78 +1134,186 @@ void SurfaceManager::createSurfaceFromPhaseMap(cv::Mat phase, CircleOutline outs
 }
 
 wavefront * SurfaceManager::readWaveFront(const QString &fileName){
-    std::ifstream file(fileName.toStdString().c_str());
-    if (!file) {
-        QString b = "Can not read file " + fileName + " " +strerror(errno);
-        QMessageBox::warning(NULL, tr("Read Wavefront File"),b);
-        return 0;
-    }
-    spdlog::get("logger")->trace("readWaveFront() step 1");
-    wavefront *wf = new wavefront();
-	wf->m_origin = WavefrontOrigin::File;
-    double width;
-    double height;
-    file >> width;
-    file >> height;
-    cv::Mat data(height,width, numType,0.);
-    spdlog::get("logger")->trace("readWaveFront() width {} height {}", width, height);
-
-    for( size_t y = 0; y < height; y++ ) {
-        for( size_t x = 0; x < width; x++ ) {
-            file >> data.at<double>(height - y-1,x);
-            //data.at<double>(height - y - 1, x) += dist(generator);
-        }
-    }
-    spdlog::get("logger")->trace("readWaveFront() step 2");
-
-    std::string line;
-    QString l;
     mirrorDlg *md = mirrorDlg::get_Instance();
-
-    double xm = (width-1)/2.,ym = (height-1)/2.,
-            radm = cv::min(xm,ym)-2 ,
-            roc = md->roc,
+    double xm,ym,radm;
+    double  roc = md->roc,
             lambda = md->lambda,
             diam = md->diameter;
-    double xo = width/2., yo = height/2., rado = 0;
-
-    std::string dummy;
-    while (getline(file, line)) {
-        l = QString::fromStdString(line);
-        std::istringstream iss(line);
-        if (l.startsWith("outside")) {
-            QStringList sl = l.split(" ");
-            xm = sl[2].toDouble();
-            radm = sl[4].toDouble();
-            ym = sl[3].toDouble();
-            continue;
-        }
-        if (l.startsWith("DIAM")){
-            iss >> dummy >> diam;
-            continue;
-        }
-        if (l.startsWith("ROC")){
-            iss >> dummy >> roc;
-            continue;
-        }
-        if (l.startsWith("Lambda")){
-            iss >> dummy >> lambda;
-            continue;
+
+    double xo, yo, rado;
+    wavefront *wf = new wavefront();
+    wf->m_origin = WavefrontOrigin::File;
+    rado=0;
+
+    if (fileName.endsWith(".npz",Qt::CaseInsensitive)){
+        //
+        // code to read npz file (WavefrontPro file)
+        //
+        double reference_wavelength=550; // wavefrontPro scales and stores the wavefront to a 550nm wavelength, not laser wavelengths (like DFTF does)
+        double obsc=0;
+        bool bAlreadyNulled=false;
+        cnpy::npz_t npz_data = cnpy::npz_load(fileName.toStdString());
+        spdlog::get("logger")->info("npz file contents");
+        bool bWavefrontLoaded = false;
+        for (const auto& element : npz_data) {
+            cnpy::NpyArray e = element.second;
+            if (element.first == "null" && e.word_size != 0)
+                bAlreadyNulled=true;  //wfpro already nulled this wavefront
+            if (e.shape.size() == 0 && e.num_vals == 1 && (e.word_size==8 || e.word_size==4)) {
+                double temp;
+                double * dval = &temp; // dval points to temp if it's an integer (32bit number)
+
+                if (e.word_size==4)
+                    temp = *(e.data<int32_t>()); // the number is a 32bit integer
+                else
+                    dval = e.data<double>(); // the number is a 64bit double
+
+                spdlog::get("logger")->info("{} size {}  word size {} num_vals {} val: {}", element.first, e.shape.size(), e.word_size, e.num_vals, *dval);
+                if (element.first == "dia")
+                    diam = *dval;
+                else if (element.first == "roc")
+                    roc = *dval;
+                //else if (element.first == "conic")
+                //    double conic = *dval;   conic is ignored so no need to read it from the npz file
+                else if (element.first == "ref_wvl")
+                    reference_wavelength= *dval;
+                else if (element.first == "laser_wvl")
+                    lambda = *dval;
+                else if (element.first == "obsc")
+                    obsc = *dval;
+                else if (element.first == "null" && std::isnan(*dval) == false)
+                    bAlreadyNulled = true;  //wfpro already nulled this wavefront
+
+            }
+            else if (e.shape.size() == 0 && e.num_vals == 1 && e.word_size==1) {
+                unsigned char * ucval = e.data<unsigned char>();
+                spdlog::get("logger")->info("{} size {}  word size {} num_vals {} val: {}", element.first, e.shape.size(), e.word_size, e.num_vals, *ucval);
+            }
+            else {
+                spdlog::get("logger")->info("{} size {}  word size {} num_vals {}", element.first, e.shape.size(), e.word_size, e.num_vals);
+                if (element.first == "wf") {
+                    if (e.shape.size() != 2 || e.word_size != 8) {
+                        delete wf;
+                        QString b = "Can not read file " + fileName + ". Wavefront not as expected (expected 2 dimensional array of doubles)";
+                        QMessageBox::warning(NULL, tr("Read Wavefront File"),b);
+                        return nullptr; // error - was expecting 2 dimensional array of doubles
+                    }
+
+                    int width = e.shape[0];
+                    int height = e.shape[1];
+                    cv::Mat data(height,width, numType,0.);
+                    for (int y=0;y<height;y++) {
+                        for (int x=0;x<width;x++) {
+                            double temp = e.data<double>()[y*width+x];
+                            if (std::isnan(temp))
+                                temp=0;
+                            data.at<double>(height-y-1,x) = temp;
+                        }
+                    }
+                    wf->data = data;
+                    if (bAlreadyNulled)
+                        wf->useSANull = false; // wavefront pro already nulled the data
+                    bWavefrontLoaded = true;
+                    radm=width/2.0 - 0.5;
+                    xm=radm;
+                    ym=radm;
+                    xo=radm;
+                    yo=radm;
+                    if (obsc!=0) {
+                        rado=obsc; // TODO: osbc may be in mm?  Or it may be in pixels.  It may be a radius or it may be a diameter.
+                    }
+
+
+                }
+            }
         }
-        if (l.startsWith("obstruction")){
-            iss >> dummy >> dummy >> xo >> yo >> rado;
-            continue;
+        if (bWavefrontLoaded == false) {
+            delete wf;
+            QString b = "Can not read file " + fileName + ". Wavefront not found";
+            QMessageBox::warning(NULL, tr("Read Wavefront File"),b);
+            return nullptr; // error - no wavefront found in npz file
         }
-        if (l.startsWith("ellipse_vertical_axis")){
-            md->m_outlineShape = ELLIPSE;
-            iss >> dummy >> md->m_verticalAxis;
+        //
+        // Regarding this following line of code. wfpro stores the data with respect to the reference wavelength (typically 550nm).  But
+        // DFTF stores the data with respect to the laser wavelength. So we need to scale the data so DFTF is happy.  Both of these values
+        // (reference_wavelength, lambda) whould always be in the npz file although an early version of wfpro didn't store the laser wavlenth
+        // there so we use the laser wavelength from the mirror dialog as a default
+        //
+        wf->data = wf->data * (reference_wavelength / lambda);
+
+    }
+    else {
+        std::ifstream file(fileName.toStdString().c_str());
+        if (!file) {
+            QString b = "Can not read file " + fileName + " " +strerror(errno);
+            QMessageBox::warning(NULL, tr("Read Wavefront File"),b);
+            delete wf;
+            return 0;
+        }
+        spdlog::get("logger")->trace("readWaveFront() step 1");
+        double width;
+        double height;
+        file >> width;
+        file >> height;
+        cv::Mat data(height,width, numType,0.);
+        spdlog::get("logger")->trace("readWaveFront() width {} height {}", width, height);
+
+        for( size_t y = 0; y < height; y++ ) {
+            for( size_t x = 0; x < width; x++ ) {
+                file >> data.at<double>(height - y-1,x);
+                //data.at<double>(height - y - 1, x) += dist(generator);
+            }
         }
-        if (l.startsWith("Do Not use null") || l.startsWith("nulled") ){
-            wf->useSANull = false;
+        spdlog::get("logger")->trace("readWaveFront() step 2");
+
+        std::string line;
+        QString l;
+
+        xm = (width-1)/2.;
+        ym = (height-1)/2.,
+        radm = cv::min(xm,ym)-2;
+        xo = width/2.;
+        yo = height/2.;
+        rado = 0;
+
+        std::string dummy;
+        while (getline(file, line)) {
+            l = QString::fromStdString(line);
+            std::istringstream iss(line);
+            if (l.startsWith("outside")) {
+                QStringList sl = l.split(" ");
+                xm = sl[2].toDouble();
+                radm = sl[4].toDouble();
+                ym = sl[3].toDouble();
+                continue;
+            }
+            if (l.startsWith("DIAM")){
+                iss >> dummy >> diam;
+                continue;
+            }
+            if (l.startsWith("ROC")){
+                iss >> dummy >> roc;
+                continue;
+            }
+            if (l.startsWith("Lambda")){
+                iss >> dummy >> lambda;
+                continue;
+            }
+            if (l.startsWith("obstruction")){
+                iss >> dummy >> dummy >> xo >> yo >> rado;
+                continue;
+            }
+            if (l.startsWith("ellipse_vertical_axis")){
+                md->m_outlineShape = ELLIPSE;
+                iss >> dummy >> md->m_verticalAxis;
+            }
+            if (l.startsWith("Do Not use null") || l.startsWith("nulled") ){
+                wf->useSANull = false;
+            }
         }
+        wf->data= data;
     }
 
-
     wf->m_outside = CircleOutline(QPointF(xm,ym), radm);
     if (rado == 0){
         xo = xm;
@@ -1292,7 +1401,6 @@ wavefront * SurfaceManager::readWaveFront(const QString &fileName){
         }
         if (rocResp == YES || messageResult == QMessageBox::Yes){
             emit rocChanged(roc);
-
         }
         else {
             roc = md->roc;
@@ -1300,7 +1408,6 @@ wavefront * SurfaceManager::readWaveFront(const QString &fileName){
         }
     }
     wf->diameter = diam;
-    wf->data= data;
     wf->roc = roc;
     wf->lambda = lambda;
     wf->wasSmoothed = false;
@@ -1329,6 +1436,7 @@ bool SurfaceManager::loadWavefront(const QString &fileName){
 
     emit enableControls(false);
     bool mirrorParamsChanged = false;
+
     std::ifstream file(fileName.toStdString().c_str());
     if (!file) {
         QString b = "Can not read file " + fileName + " " +strerror(errno);