Index: _PolConvert.cpp
===================================================================
--- _PolConvert.cpp	(revision 9901)
+++ _PolConvert.cpp	(working copy)
@@ -1,5 +1,4 @@
 /*
-
 # Copyright (c) Ivan Marti-Vidal 2015-2020. 
 #               EU ALMA Regional Center. Nordic node.
 #               Observatori Astronomic, Universitat de Valencia
@@ -52,7 +51,12 @@
 
 */
  
+//z old code 1, new code 0
+#define IVAN_OLD    1
 
+//z straight up merge
+#define IVAN_MERGE  0
+
 #include <Python.h>
 
 // compiler warning that we use a deprecated NumPy API
@@ -316,7 +320,8 @@
   double *AntCoordArr = (double *)PyArray_DATA(antcoordObj);
   int *AntMountArr = (int *)PyArray_DATA(antmountObj);
   double *SouCoordArr = (double *)PyArray_DATA(PyList_GetItem(soucoordObj,1));
-    
+//z
+  double *SouCoordRA = (double *)PyArray_DATA(PyList_GetItem(soucoordObj,0));
   Geometry->NtotSou = (int) PyArray_DIM(PyList_GetItem(soucoordObj,1),0);
   Geometry->NtotAnt = (int) PyArray_DIM(antcoordObj,0);
 
@@ -328,6 +333,8 @@
   Geometry->BaseLine[2] = new double[Nbas+Geometry->NtotAnt+1];
   Geometry->SinDec = new double[Geometry->NtotSou];
   Geometry->CosDec = new double[Geometry->NtotSou];
+//z:
+  Geometry->RA = new double[Geometry->NtotSou];
   Geometry->AntLon = new double[Geometry->NtotAnt];
   Geometry->Mount = new int[Geometry->NtotAnt];
   Geometry->Lat = new double[Geometry->NtotAnt];
@@ -364,8 +371,11 @@
   };
 
   for (i=0; i<Geometry->NtotSou;i++){
-      Geometry->SinDec[i] = sin(SouCoordArr[i+Geometry->NtotSou]);
-      Geometry->CosDec[i] = cos(SouCoordArr[i+Geometry->NtotSou]);
+//z      Geometry->SinDec[i] = sin(SouCoordArr[i+Geometry->NtotSou]);
+//z      Geometry->CosDec[i] = cos(SouCoordArr[i+Geometry->NtotSou]);
+        Geometry->SinDec[i] = sin(SouCoordArr[i]);
+        Geometry->CosDec[i] = cos(SouCoordArr[i]);
+        Geometry->RA[i] = SouCoordRA[i];
 //      std::cout<< i << " "<< SouCoordArr[i]*180./3.1415926535 << "\n";
   };
 
@@ -587,13 +597,32 @@
 
 // Which IFs do we convert?
   int IFs2Conv[nIFconv];
+//z
+#if IVAN_OLD
+#warning "IVAN_OLD CODE selected A"
   for (ii=0; ii<nIFconv; ii++) {
-      if (doAll){IFs2Conv[ii]=ii+1;} else {
-        IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii));
+    if (doAll){IFs2Conv[ii]=ii+1;} else {
+      IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii)); };
+  }
+#else // IVAN_OLD
+#warning "IVAN_NEW CODE selected A"
+  // looks like a shift by 1 here
+  for (ii=0; ii<nIFconv; ii++) {
+    if (doAll){IFs2Conv[ii]=ii;} else {
+      IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii)) - 1; };
+  }
+#endif // IVAN_OLD
+#if IVAN_MERGE
+#warning "IVAN_MERGED CODE selected A"
+  for (ii=0; ii<nIFconv; ii++) {
+//z      if (doAll){IFs2Conv[ii]=ii+1;} else {
+//z        IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii));
+      if (doAll){IFs2Conv[ii]=ii;} else {
+        IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii)) - 1;
       };
   }; 
+#endif // IVAN_MERGE
 
-
   int nnu = DifXData->getNfreqs();
   int ALMARefAnt = -1;
   // If no calAPP is used, do not look for any extra X-Y phase offset.
@@ -719,6 +748,9 @@
 
 // Prepare plotting files:
   int noI = -1;
+//z
+#if IVAN_OLD
+#warning "IVAN_OLD CODE selected B"
   for (ii=0; ii<nIFplot; ii++){
       sprintf(message,"POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i",IFs2Plot[ii]+1);
       printf("Writing %s\n", message);
@@ -729,6 +761,36 @@
          fwrite(&noI,sizeof(int),1,plotFile[ii]);
       };
   };
+#else // IVAN_OLD
+#warning "IVAN_NEW CODE selected B"
+  for (ii=0; ii<nIFconv; ii++) {
+      sprintf(message,"POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i",IFs2Conv[ii]+1);
+      printf("Writing %s\n", message);
+      plotFile[ii] = fopen(message,"wb");
+      if (IFs2Conv[ii]>=0 && IFs2Conv[ii]<nnu){
+         fwrite(&nchans[IFs2Conv[ii]],sizeof(int),1,plotFile[ii]);
+      } else {
+         fwrite(&noI,sizeof(int),1,plotFile[ii]);
+      };
+  };
+#endif // IVAN_OLD
+#if IVAN_MERGE
+#warning "IVAN_MERGED CODE selected B"
+//z  for (ii=0; ii<nIFplot; ii++){
+//z      sprintf(message,"POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i",IFs2Plot[ii]+1);
+  for (ii=0; ii<nIFconv; ii++) {
+      sprintf(message,"POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i",IFs2Conv[ii]+1);
+      printf("Writing %s\n", message);
+      plotFile[ii] = fopen(message,"wb");
+//z      if (IFs2Plot[ii]>=0 && IFs2Plot[ii]<nnu){
+//z         fwrite(&nchans[IFs2Plot[ii]],sizeof(int),1,plotFile[ii]);
+      if (IFs2Conv[ii]>=0 && IFs2Conv[ii]<nnu){
+         fwrite(&nchans[IFs2Conv[ii]],sizeof(int),1,plotFile[ii]);
+      } else {
+         fwrite(&noI,sizeof(int),1,plotFile[ii]);
+      };
+  };
+#endif // IVAN_MERGE
  
 
   if(doNorm){
@@ -758,7 +820,20 @@
 
   for (currAntIdx=0; currAntIdx<nALMA; currAntIdx++) {
       for (im=0; im<nIFconv; im++) {
+//z
+#if IVAN_OLD
+#warning "IVAN_OLD CODE selected C"
         ii = IFs2Conv[im] - 1;
+#else // IVAN_OLD
+#warning "IVAN_NEW CODE selected C"
+        ii = IFs2Conv[im];
+#endif // IVAN_OLD
+#if IVAN_MERGE
+#warning "IVAN_MERGED CODE selected C"
+//z     ii = IFs2Conv[im] - 1;
+        ii = IFs2Conv[im];
+#endif // IVAN_MERGE
+//z
         for (ij=0; ij<nchans[ii]; ij++){
           PrioriGains[currAntIdx][im][ij] *=
             DifXData->getAmpRatio(currAntIdx, im, ij);
@@ -780,13 +855,49 @@
 //  for (ii=0; ii<nnu; ii++) open-brace
   for (im=0; im<nIFconv; im++) {
 
+//z
+#if IVAN_OLD
+#warning "IVAN_OLD CODE selected D"
     ii = IFs2Conv[im] - 1;
+#else // IVAN_OLD
+#warning "IVAN_NEW CODE selected D"
+    ii = IFs2Conv[im];
+#endif // IVAN_OLD
+#if IVAN_MERGE
+#warning "IVAN_MERGED CODE selected D"
+//z ii = IFs2Conv[im] - 1;
+    ii = IFs2Conv[im];
+#endif // IVAN_MERGE
+//z
 
+//z
+#if IVAN_OLD
+#warning "IVAN_OLD CODE selected E"
     bool plotIF = false;
     int IFplot = 0;
     for (ij=0; ij<nIFplot; ij++){
-      if (IFs2Plot[ij]==ii){plotIF=true;IFplot=ij; break;};
+      if (IFs2Plot[ij]==ii){
+        plotIF=true;
+        IFplot=ij; break;};
     };
+#else // IVAN_OLD
+#warning "IVAN_NEW CODE selected E"
+    int IFplot = 0;
+    for (ij=0; ij<nIFplot; ij++){
+      if (IFs2Plot[ij]==ii){
+        IFplot=ij; break;};
+#endif // IVAN_OLD
+#if IVAN_MERGE
+#warning "IVAN_MERGED CODE selected E"
+//zz    bool plotIF = false;
+    int IFplot = 0;
+    for (ij=0; ij<nIFplot; ij++){
+      if (IFs2Plot[ij]==ii){
+//zz    plotIF=true;
+        IFplot=ij; break;};
+    };
+#endif // IVAN_MERGE
+//zz
 
  //   if (ii >= nnu) {
  //     sprintf(message,"ERROR! DATA DO NOT HAVE IF #%i !!\n",ii);  
@@ -1246,10 +1357,32 @@
 
 // Shall we write in plot file?  auxB -> auxB2 to avoid confusion
         //indent level within time range
+
+//z
+//z QUESTION: so why does plotIF disappear here?
+//z note auxB -> auxB2
+#if IVAN_OLD
+#warning "IVAN_OLD CODE selected F"
         auxB2 = (currT>=plRange[0] && currT<=plRange[1] &&
                 plotIF && (calField<0 || currF==calField));
                 //plAnt == otherAnt);
+#else // IVAN_OLD
+#warning "IVAN_NEW CODE selected F"
+        auxB2 = (currT>=plRange[0] && currT<=plRange[1] &&
+                (calField<0 || currF==calField));
+                //plAnt == otherAnt);
+#endif // IVAN_OLD
+#if IVAN_MERGE
+#warning "IVAN_MERGED CODE selected F"
+//z     auxB2 = (currT>=plRange[0] && currT<=plRange[1] &&
+//z             plotIF && (calField<0 || currF==calField));
+                //plAnt == otherAnt);
+        auxB2 = (currT>=plRange[0] && currT<=plRange[1] &&
+                (calField<0 || currF==calField));
+                //plAnt == otherAnt);
+#endif // IVAN_MERGE
 
+
         //indent level within time range
 // Convert:
         if(Phased){
Index: _getAntInfo.cpp
===================================================================
--- _getAntInfo.cpp	(revision 9901)
+++ _getAntInfo.cpp	(working copy)
@@ -72,17 +72,25 @@
     "Returns antenna coordinates";
 static char getMounts_docstring[] =
     "Returns the mount types";
+//z
+//static char getNames_docstring[] = 
+//    "Returns the antenna names (codes)";
 
 /* Available functions */
 static PyObject *getAntInfo(PyObject *self, PyObject *args);
 static PyObject *getCoords(PyObject *self, PyObject *args);
 static PyObject *getMounts(PyObject *self, PyObject *args);
+//z
+//static PyObject *getNames(PyObject *self, PyObject *args);
 
+
 /* Module specification */
 static PyMethodDef module_methods[] = {
     {"getAntInfo", getAntInfo, METH_VARARGS, getAntInfo_docstring},
     {"getCoords", getCoords, METH_VARARGS, getCoords_docstring},
-    {"getMounts", getMounts, METH_VARARGS, getMounts_docstring},
+//z {"getMounts", getMounts, METH_VARARGS, getMounts_docstring},
+    {"getMounts", getMounts, METH_VARARGS, getMounts_docstring}, //,
+//  {"getNames", getNames, METH_VARARGS, getNames_docstring},
     {NULL, NULL, 0, NULL}   /* terminated by list of NULLs, apparently */
 };
 
@@ -216,7 +224,8 @@
   fits_close_file(ifile, &status);
 
   if(status){
-    printf("\n\nPROBLEM LOSING FITS-IDI!  ERR: %i\n\n",status);
+//z printf("\n\nPROBLEM LOSING FITS-IDI!  ERR: %i\n\n",status);
+    printf("\n\nPROBLEM CLOSING FITS-IDI!  ERR: %i\n\n",status);
     return Py_BuildValue("i",3);
   };
 
@@ -226,33 +235,60 @@
 };
 
 
+
+
+
+
+
 static PyObject *getCoords(PyObject *self, PyObject *args){
 
 // Build numpy array:
 
 PyObject *CoordArr;
-long CD[2] = {Nants,3};
+//z long CD[2] = {Nants,3};
+//z
+//z GBC had already removed Py_INCREF.
+//z Py_INCREF(CoordArr); -- original code
+//z CoordArr = PyArray_SimpleNewFromData(2, &CD[0], NPY_DOUBLE, (void*) Coords);
+//z Py_INCREF(CoordArr); -- needed only if we retain it and use it
+int nd = 2;
+npy_intp* CD = new npy_intp[nd];
 
-//Py_INCREF(CoordArr); -- original code
-CoordArr = PyArray_SimpleNewFromData(2, &CD[0], NPY_DOUBLE, (void*) Coords);
-//Py_INCREF(CoordArr); -- needed only if we retain it and use it
+CD[0] = Nants; CD[1] = 3; // = {Nants,3};
 
+//Py_INCREF(CoordArr);
+CoordArr = PyArray_SimpleNewFromData(nd, CD, NPY_DOUBLE, (void*) Coords);
+
 return CoordArr;
 
 };
 
 
+
+
+
+
+
 static PyObject *getMounts(PyObject *self, PyObject *args){
 
+
 // Build numpy array:
 PyObject *MountArr;
 
-long MD[1] = {Nants};
+//z long MD[1] = {Nants};
+int nd = 1;
+npy_intp* MD = new npy_intp[nd];
+MD[0] = Nants;
 
-//Py_INCREF(MountArr); -- original code
-MountArr = PyArray_SimpleNewFromData(1, &MD[0], NPY_INT, (void*) Mounts);
-//Py_INCREF(MountArr); -- needed only if we retain it and use it
+//z GBC had already removed Py_INCREF
+//z Py_INCREF(MountArr); -- original code
+//z MountArr = PyArray_SimpleNewFromData(1, &MD[0], NPY_INT, (void*) Mounts);
+//z Py_INCREF(MountArr); -- needed only if we retain it and use it
 
+//z
+//Py_INCREF(MountArr);
+MountArr = PyArray_SimpleNewFromData(nd, MD, NPY_INT, (void*) Mounts);
+
 return MountArr;
 
 };
Index: DataIO.cpp
===================================================================
--- DataIO.cpp	(revision 9901)
+++ DataIO.cpp	(working copy)
@@ -106,42 +106,63 @@
 
 
 
-void DataIO::getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &P1, double &P2){
+//z void DataIO::getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &P1, double &P2){
+void DataIO::getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &MJD, double &P1, double &P2){
 
-double V2, Bx, By; //, Bz;
-double CH, SH, CT1, CT2, HAng, H1, H2, Elev1, Elev2;
-int ibas;
+//z double V2, Bx, By; //, Bz;
+//z double CH, SH, CT1, CT2, HAng, H1, H2, Elev1, Elev2;
+//zz double V2, Bx, By, GMST; //, Bz;
+double GMST;
+//double CDec, SDec, SRA, TLat1,TLat2, Lon1,Lon2,CH, SH;
+double CT1, CT2, HAng, H1, H2, Elev1, Elev2;
+//zz int ibas;
 
 Elev1 = 0.0; Elev2 = 0.0;
 
+//z:
+double days = MJD/86400.;
+double t = (days-51544.0)/36525.;
+double Hh = days - floor(days);
+double GMsec = 24110.54841 + 8640184.812866*t + 0.093104*t*t - 0.0000062*t*t*t;
+GMST = (GMsec/86400. + Hh)*2.*3.1415926535;
+//z:
+//CDec = Geometry->SinDec[sidx];
+//SDec = Geometry->CosDec[sidx];
+//TLat1 = tan(Geometry->Lat[Ant1]);
+//TLat2 = tan(Geometry->Lat[Ant1]);
+//Lon1 = Geometry->AntLon[Ant1];
+//Lon2 = Geometry->AntLon[Ant2];
+//SRA = Geometry->RA[sidx];
 
 // Dummy value if autocorrelation:
 if (Ant1==Ant2){P1 = -1.e9; P2 = -1.e9; return;};
 
+//z more blank lines
 
 if(sidx<Geometry->NtotSou && Ant1<Geometry->NtotAnt && Ant2<Geometry->NtotAnt){
 
-  V2 = Geometry->SinDec[sidx]*UVW[1] - Geometry->CosDec[sidx]*UVW[2];
-  ibas = Geometry->BasNum[Ant1][Ant2];
+//z  V2 = Geometry->SinDec[sidx]*UVW[1] - Geometry->CosDec[sidx]*UVW[2];
+//zz   ibas = Geometry->BasNum[Ant1][Ant2];
 
+//z  if (ibas<0){
+//z  ibas = -ibas;
+//z   Bx = -Geometry->BaseLine[0][ibas];
+//z   By = -Geometry->BaseLine[1][ibas];
+//z// Bz = -Geometry->BaseLine[2][ibas];
+//z  } else {
+//z   Bx = Geometry->BaseLine[0][ibas];
+//z   By = Geometry->BaseLine[1][ibas];
+//z// Bz = Geometry->BaseLine[2][ibas];
+//z  };
+//z
+//z  CH = (UVW[0]*By - V2*Bx); // /(By**2. + Bx**2.);
+//z  SH = (UVW[0]*Bx + V2*By); // /(By**2. + Bx**2.);
 
-  if (ibas<0){
-  ibas = -ibas;
-   Bx = -Geometry->BaseLine[0][ibas];
-   By = -Geometry->BaseLine[1][ibas];
-//   Bz = -Geometry->BaseLine[2][ibas];
-  } else {
-   Bx = Geometry->BaseLine[0][ibas];
-   By = Geometry->BaseLine[1][ibas];
-//   Bz = Geometry->BaseLine[2][ibas];
-  };
-
-  CH = (UVW[0]*By - V2*Bx); // /(By**2. + Bx**2.);
-  SH = (UVW[0]*Bx + V2*By); // /(By**2. + Bx**2.);
   CT1 = Geometry->CosDec[sidx]*tan(Geometry->Lat[Ant1]);
   CT2 = Geometry->CosDec[sidx]*tan(Geometry->Lat[Ant2]);
 
-  HAng = atan2(SH,CH);
+//z  HAng = atan2(SH,CH);
+  HAng = GMST - Geometry->RA[sidx];
   H1 = HAng + Geometry->AntLon[Ant1];
   H2 = HAng + Geometry->AntLon[Ant2];
 
@@ -148,10 +169,11 @@
   // sin(lat)*sin(dec)+cos(lat)cos(dec)cos(H)
 
   if (Geometry->Mount[Ant1] > 3){
-  Elev1 = asin(sin(Geometry->Lat[Ant1])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant1])*Geometry->CosDec[sidx]*cos(H1));};
+   Elev1 = asin(sin(Geometry->Lat[Ant1])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant1])*Geometry->CosDec[sidx]*cos(H1));};
 
   if (Geometry->Mount[Ant2] > 3){
-  Elev2 = asin(sin(Geometry->Lat[Ant1])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant1])*Geometry->CosDec[sidx]*cos(H1));};
+//z  Elev2 = asin(sin(Geometry->Lat[Ant1])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant1])*Geometry->CosDec[sidx]*cos(H1));};
+   Elev2 = asin(sin(Geometry->Lat[Ant2])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant2])*Geometry->CosDec[sidx]*cos(H2));};
  
 
   switch (Geometry->Mount[Ant1]){
@@ -178,6 +200,8 @@
 
 //  P2 = atan2(sin(H2), CT2 - Geometry->SinDec[sidx]*cos(H2));
 
+//z:
+//if(Ant1==1){printf("%i  %i  %.16e  %.3f  %.3f  %.3f  %.3f\n",Ant1, Ant2, MJD, GMST, HAng, P1, P2);};
 
 } else {
 
Index: _PolGainSolve.cpp
===================================================================
--- _PolGainSolve.cpp	(revision 9901)
+++ _PolGainSolve.cpp	(working copy)
@@ -1654,6 +1654,7 @@
 
        if (af1>=0){
          RateResVec[af1] += Weights[i][BNum]*(BLRates00[i][BNum] + BLRates11[i][BNum])/2.;
+//z      // printf("%i | %i-%i %.2e  %.2e \n",i, CalAnts[a1],CalAnts[a2], BLRates00[i][BNum],BLRates11[i][BNum]);
          DelResVec00[af1] += Weights[i][BNum]*BLDelays00[i][BNum];
          DelResVec11[af1] += Weights[i][BNum]*BLDelays11[i][BNum];
          Hessian[af1*NantFit + af1] += Weights[i][BNum];
@@ -1701,21 +1702,36 @@
 
 
 printf("\n\nHessian Globalization Matrix:\n\n");
-bool isSingular, tempSing;
+//zz bool isSingular, tempSing;
+bool isSingular;
 isSingular=false;
 for (i=0; i<NantFit; i++){
   printf("  ");
-  tempSing = true;
+//z  tempSing = true;
+  if (Hessian[i*NantFit+i]==0.0){isSingular=true;};
   for (j=0; j<NantFit; j++){
-     if (Hessian[i*NantFit+j]!=0.0){tempSing=false;};
+//z     if (Hessian[i*NantFit+j]!=0.0){tempSing=false;};
      printf("%+.2e ",Hessian[i*NantFit+j]);
   };
-  if (tempSing){isSingular=true;};
+//z  if (tempSing){isSingular=true;};
+//z printf("\n");
+//z };
 printf("\n");
 };
+if(isSingular){printf("Possible singular matrix!\n");};
+printf("\n");
+
+printf("\n\n Residual baseline phase quantities:\n\n");
+for (i=0; i<NantFit; i++){
+  printf("  ");
+  printf(" %.2e | %.2e | %.2e\n",DelResVec00[i],DelResVec11[i],RateResVec[i]);
+};
+
+// GBC added these which may now be redundant
 if (!isSingular) printf("\n(not Singular)\n");
 else             printf("\n(is  Singular)\n");
 
+printf("\n");
 
 // The Hessian's inverse can be reused for rates, delays and phases!
 gsl_death_by = GSL_SUCCESS;
@@ -1753,6 +1769,8 @@
 }
 
 if(!isSingular){
+//z     printf("Globalizing solutions\n");
+        printf("Globalizing solutions: ");
         printf("doing gsl_linalg_LU_solve 3x...\n");
 	gsl_linalg_LU_solve (&mm.matrix, permm, &RateInd.vector, xx);
 	gsl_linalg_LU_solve (&mm.matrix, permm, &Del00Ind.vector, dd0);
@@ -2274,6 +2292,8 @@
 
   for(k=0;k<NBas;k++){Tm[k]=Times[currIF][0];};
 
+//z
+// printf("NVis: %i\n",NVis[currIF]); fflush(stdout);
 
   for (k=0; k<NVis[currIF]; k++){
 
@@ -2436,6 +2456,9 @@
 
     if(BNum>=0){
 
+//z
+    //  printf("Fitting for %i %i\n",a1,a2); fflush(stdout);
+
     AvVis[BNum] += 1;
 
     FeedFactor1 = std::polar(1.0, feedAngle[a1-1])*PA1[currIF][k]; 
Index: task_polconvert.py
===================================================================
--- task_polconvert.py	(revision 9901)
+++ task_polconvert.py	(working copy)
@@ -51,7 +51,7 @@
 from __future__ import absolute_import
 from __future__ import print_function
 __version__ = "1.8.1"
-date = 'Oct 11 2020'     
+date = 'Feb 09 2021'     
 
 
 ################
@@ -112,15 +112,24 @@
 import pickle as pk
 
 if sys.version_info.major < 3:
+  try:
     from taskinit import *
     ms = gentools(['ms'])[0]
     tb = gentools(['tb'])[0]
+  except Exception as ex:
+    print('unable to load casa tools in python2\n\n')
+    raise ex
 else:
+  try:
     # from taskutil import *
     from casatools import ms as ms_casa
     from casatools import table as tb_casa
     ms = ms_casa()
     tb = tb_casa()
+  except Exception as ex:
+    print('unable to load casa tools in python3\n\n')
+    raise ex
+    
 
 #########################################################
 ###########################
@@ -134,25 +143,41 @@
 
 if __name__=='__main__':
 
+##z not changed yet
+  print('setting up variables')
+
+  # grep name= polconvert.xml | cut -d\" -f 4 | tr \\012 ,
+  global IDI,OUTPUTIDI,DiFXinput,DiFXcalc,doIF,linAntIdx,Range,ALMAant
+  global spw,calAPP,calAPPTime,APPrefant,gains,interpolation,gainmode
+  global XYavgTime,dterms,amp_norm,XYadd,XYdel,XYratio,usePcal,swapXY
+  global swapRL,feedRotation,IDI_conjugated,plotIF,plotRange,plotAnt
+  global excludeAnts,excludeBaselines,doSolve,solint,doTest,npix,solveAmp
+  global solveMethod,calstokes,calfield
+
   # set things that the task machinery should set 
   taskname           = "polconvert"
-  IDI                =  "bm494e-0-b1_1200.difx"
-  OUTPUTIDI          =  "1200_POLCONVERT.difx"
-  DiFXinput          =  "bm494e-0-b1_1200.input"
-  DiFXcalc           =  "bm494e-0-b1_1200.calc"
-  doIF               =  [35, 36, 37, 38]
+  IDI                =  "no-such-dir.difx"
+  OUTPUTIDI          =  "no-such-dir.difx"
+  DiFXinput          =  "no-such-file.input"
+  DiFXcalc           =  "no-such-file.calc"
+  doIF               =  [6660,6661]
   linAntIdx          =  [1]
   Range              =  []
-  ALMAant            =  "APP_DERIVERABLES/SPECLINE_0.concatenated.ms.ANTENNA"
+  ALMAant            =  "no-such-dir/label.concatenated.ms.ANTENNA"
   spw                =  0
-  calAPP             =  "APP_DERIVERABLES/SPECLINE_0.concatenated.ms.calappphase"
+  calAPP             =  "no-such-dir/label.concatenated.ms.calappphase"
   calAPPTime         =  [0.0, 5.0]
   APPrefant          =  ""
-  gains              =  [['APP_DERIVERABLES/SPECLINE_0.concatenated.ms.bandpass-zphs', 'APP_DERIVERABLES/SPECLINE_0.concatenated.ms.flux_inf.APP', 'APP_DERIVERABLES/SPECLINE_0.concatenated.ms.phase_int.APP.XYsmooth', 'APP_DERIVERABLES/SPECLINE_0.calibrated.ms.XY0.APP', 'APP_DERIVERABLES/SPECLINE_0.calibrated.ms.Gxyamp.ALMA']]
+  gains              =  [[
+    'no-such-dir/label.concatenated.ms.bandpass-zphs',
+    'no-such-dir/label.concatenated.ms.flux_inf.APP',
+    'no-such-dir/label.concatenated.ms.phase_int.APP.XYsmooth',
+    'no-such-dir/label.calibrated.ms.XY0.APP',
+    'no-such-dir/label.calibrated.ms.Gxyamp.ALMA']]
   interpolation      =  []
   gainmode           =  [['G', 'T', 'G', 'G', 'G']]
   XYavgTime          =  0.0
-  dterms             =  ['APP_DERIVERABLES/SPECLINE_0.calibrated.ms.Df0gen.ALMA']
+  dterms             =  ['no-such-dir/label.calibrated.ms.Df0gen.ALMA']
   amp_norm           =  0.03
   XYadd              =  {}
   XYdel              =  {}
@@ -166,6 +191,7 @@
   plotRange          =  [0, 0, 0, 0, 2, 0, 0, 0]
   plotAnt            =  2
   excludeAnts        =  []
+  excludeBaselines   =  []
   doSolve            =  -1
   solint             =  [1, 1]
   doTest             =  True
@@ -174,6 +200,22 @@
   solveMethod        =  "gradient"
   calstokes          =  [1.0, 0.0, 0.0, 0.0]
   calfield           =  -1
+
+  try:
+    # rename polconvert.last to polconvert.test to use it
+    pctest = open('polconvert.test')
+    changes = pctest.readlines()
+    pctest.close()
+    exec('\n'.join(changes))
+    print('updated defaults from polconvert.test')
+  except Exception as ex:
+    print('you cannot run this stand-alone without a polconvert.test\n' +
+          'file that re-assigns the variables.  You can use a previous\n' +
+          'polconvert.last if you copy/rename it to polconvert.test\n')
+    raise ex
+
+##z no changes above, yet
+
 #
 #
 #
@@ -714,15 +756,21 @@
     printError("ERROR! OUTPUTIDI should be a string!")
 
   if type(plotIF) is int:
-    plotIF = [plotIF]
+##z
+#   plotIF = [plotIF]
+    if plotIF >0:
+      plotIF = [plotIF]
+    else:
+      plotIF = []
   for pli in plotIF:
     if type(pli) is not int:
       printError("ERROR! plotIF should be an integer or a list of integers!")
 
-  try:
-    plotAnt = int(plotAnt)
-  except:
-    printError("ERROR! plotAnt should be an integer!")
+##z
+# try:
+#   plotAnt = int(plotAnt)
+# except:
+#   printError("ERROR! plotAnt should be an integer!")
 
 
   try:
@@ -857,10 +905,10 @@
   except:
     printError("Bad format for calAPPTime. Should be a list of 2 floats!")
 
+##z
+# if plotAnt in linAntIdx:
+#   printMsg("WARNING: Plotting will involve autocorrelations. \nThis has not been fully tested!") 
 
-  if plotAnt in linAntIdx:
-    printMsg("WARNING: Plotting will involve autocorrelations. \nThis has not been fully tested!") 
-
   printMsg("XYadd is %s"%str(XYadd))
   printMsg("XYdel is %s"%str(XYdel))
   printMsg("XYratio is %s"%str(XYratio))
@@ -1044,6 +1092,8 @@
       if decappUnit=='DEGREES':
         soucoords[1] *= np.pi/180.
 
+##z
+      antcodes = [ff[:2] for ff in ffile['ANTENNA'].data['ANNAME']]
       ffile.close()
 
 # THESE LINES FAIL IF ORBPARM IS PRESENT IN ARRAY GEOMETRY!
@@ -1057,7 +1107,7 @@
       else:
         antcoords = gA.getCoords()
         antmounts = gA.getMounts()
-        antcodes = ['%02i'%i for i in range(1,len(antmounts)+1)]
+##z #   antcodes = ['%02i'%i for i in range(1,len(antmounts)+1)]
     except:
       printMsg('WARNING! This FITS-IDI file has missing information!\nPolConvert may not calibrate properly.')
   else:
@@ -1068,6 +1118,8 @@
 ######
 
 
+
+
 ######
 # IF THIS IS A SWIN DATASET, READ THE INPUT FILE INTO 
 # A METADATA LIST:
@@ -1140,6 +1192,7 @@
  
 
   else:
+
 # READ FREQUENCY INFO TO HELP SELECTING THE SPW AUTOMATICALLY:
     import pyfits as pf
     fitsf = pf.open(IDI)
@@ -1149,13 +1202,30 @@
     IFchan = nch
     Nr = fitsf['FREQUENCY'].header['NO_BAND']
     sgn = {True:1.0,False:-1.0}[bw>0.0]
-    FrInfo = {'FREQ (MHZ)':[nu0/1.e6], 'BW (MHZ)':[bw*nch/1.e6], 'SIGN':[sgn]}
+##z
+##z FrInfo = {'FREQ (MHZ)':[nu0/1.e6], 'BW (MHZ)':[bw*nch/1.e6], 'SIGN':[sgn]}
+    FrInfo = {'FREQ (MHZ)':[], 'BW (MHZ)':[], 'SIGN':[], 'NUM CHANNELS':[]}
+    if sgn:
+      FrInfo['SIDEBAND'] = ['U' for i in range(Nr)]
+    else:
+      FrInfo['SIDEBAND'] = ['L' for i in range(Nr)]
+
+    metadata = []
+##z
     for i in range(Nr):
-      FrInfo['FREQ (MHZ)'] += [(nu0 + bw*nch/1.e6)/1.e6]
+##z   FrInfo['FREQ (MHZ)'] += [(nu0 + bw*nch/1.e6)/1.e6]
+      FrInfo['FREQ (MHZ)'] += [(nu0 + i*bw*nch)/1.e6]
       FrInfo['BW (MHZ)'] += [bw*nch/1.e6]
       FrInfo['SIGN'] += [sgn]
       FrInfo['NUM CHANNELS'] += [int(nch)]
+      freqs = nu0 + np.linspace((sgn-1.)/2.,(sgn+1.)/2.,nch,endpoint=False)*bw
+      metadata.append(freqs)
 
+    FrInfo['CHANS TO AVG'] = [1 for i in range(Nr)]
+    FrInfo['OVERSAMPLE FAC.'] = [1 for i in range(Nr)]
+    FrInfo['DECIMATION FAC.']=[1 for i in range(Nr)]
+
+
     if len(doIF)==0:
      doIF = list(range(1,1+fitsf['FREQUENCY'].header['NO_BAND']))
 
@@ -1163,6 +1233,7 @@
 
 
 
+
 # ANTENNAS TO PARTICIPATE IN THE GAIN ESTIMATES:
   nTotAnt = len(antcoords)
 
@@ -1169,16 +1240,41 @@
   calAnts = []
   for exA in antcodes:
     if exA not in excludeAnts:
+##z   #print exA, excludeAnts, exA in excludeAnts
       calAnts.append(antcodes.index(exA)+1)
       ### = [i+1 for i in range(len(antcoords)) if i+1 not in excludeAnts]
     else:
       printMsg("Excluding antenna %s"%str(exA))
 
+##z following section
+  try:
+    plotAnt = int(plotAnt)
+  except:
+    if plotAnt not in antcodes:
+      printError("Reference antenna %s is not found in metadata!"%str(plotAnt))
+    else:
+      plotAnt = antcodes.index(plotAnt)+1
 
+  for i in range(len(linAntIdx)):
+    try:
+      linAntIdx[i] = int(linAntIdx[i])
+    except:
+      if linAntIdx[i] not in antcodes:
+        linAntIdx[i] = antcodes.index(linAntIdx[i])+1
+
+
+  if plotAnt in linAntIdx:
+    printMsg(
+        "WARNING: Plotting will involve autocorrelations. \nThis has not been fully tested!")
+
+
+
+
+
   FlagBas1 = []
   FlagBas2 = []
   for fbi in excludeBaselines:
-    print(fbi, antcodes)
+##z # print(fbi, antcodes)
     printMsg("Excluding baseline %s"%str(fbi))
     if fbi[0] in antcodes and fbi[1] in antcodes:
       FlagBas1.append(antcodes.index(fbi[0])+1) ### = np.array([int(i[0]+1) for i in excludeBaselines])
@@ -1211,7 +1307,6 @@
 
 
 
-
 #######################
 ##### GET SPECTRAL WINDOW AUTOMATICALLY:
   if isPhased and spw < 0:
@@ -1401,6 +1496,7 @@
    except:
      printError("Bad time range format for plotRange!")
 
+
   if len(Range) == 0:
     Ran = np.array([0.,1.e20])
   else:
@@ -1910,7 +2006,7 @@
 
 
   if DEBUG:
-    PC_Params = [nALMATrue, doIF, plotAnt, len(allants), doIF,
+    PC_Params = [nALMATrue, plotIF, plotAnt, len(allants), doIF,
         swapXY, ngain, NSUM, kind, len(gaindata), len(dtdata), OUTPUT,
         linAntIdxTrue, plRan, Ran, allantidx, len(nphtimes), len(antimes),
         len(refants), len(asdmtimes),  doTest, doSolve, doConj, doAmpNorm,
@@ -1935,6 +2031,7 @@
 
   printMsg("\n###\n### Done with PolConvert (status %d).\n###" % (didit))
 
+##z # raw_input('HOLD')
 
   if didit != 0:
     printError("\n\n ERROR IN POLCONVERT!\n\n")
@@ -2136,8 +2233,8 @@
 
        i += 1
 
-       while Chi2_0>currChi2: 
-
+##z    while Chi2_0>currChi2: 
+       while Chi2_0>=currChi2:
          i += 1     
          LMTune *= KFacRaise
          ptst0 = np.copy(currP)
@@ -2148,8 +2245,13 @@
          if i>=MAXIT:
            break
 
-       relchange = (currChi2 - Chi2_0)/Chi2_0
+##z    relchange = (currChi2 - Chi2_0)/Chi2_0
+       if Chi2_0 >0.0:
+         relchange = (currChi2 - Chi2_0)/Chi2_0
+       else:
+         printError("\n\n  Problem in PolGainSolve.\n")         
 
+
 # No improvement:
        if currChi2<minChi2:
          minChi2 = currChi2
@@ -2217,7 +2319,10 @@
 
      Chi2_final = PS.GetChi2(minGains,LMTune,Ch0,Ch1,1)
 
-     FLIP = Chi2_final < 0.0  # Flip gains by 180 degrees.
+##z  FLIP = Chi2_final < 0.0  # Flip gains by 180 degrees.
+     FLIP = False # Chi2_final < 0.0  # Flip gains by 180 degrees.
+
+# GBC debugging:
      if FLIP: sys.stdout.write(' Flipped\n')
      else:    sys.stdout.write(' NotFlip\n')
      printMsg("    Final error: %.3e in ChSq"%(np.abs(relchange)))
@@ -2456,7 +2561,8 @@
       sub2.set_xlim((np.min(Freq2Plot) - Dnu*0.1,np.max(Freq2Plot) + Dnu*0.45))
       sub2.set_ylim((0.,2.5))
 
-      sub2.legend(numpoints=1)
+##z   sub2.legend(numpoints=1)
+      sub1.legend(numpoints=1)
       sub1.set_ylabel('Cross-Phase (deg.)')
       sub2.set_ylabel('Cross-Amp (Norm.)')
       sub2.set_xlabel('Frequency (GHz)')
@@ -2536,7 +2642,8 @@
      if os.stat("POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i"%pli).st_size>10:
        GoodIFs.append(pli)
      else:
-       printMsg("WARNING! IF %i was NOT polconverted properly\n"%pli)
+       printMsg(("WARNING! IF %i was NOT polconverted properly\n"%pli) +
+         ("POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i missing"%pli))
 
 # start of GoodIFs pli loop
    for pli in GoodIFs:
@@ -3010,6 +3117,9 @@
     printMsg(str(ex))
   ofile.close()
   printMsg('PolConvert.XYGains.dat was written with CGains' + str(CGains.keys()))
+
+# raw_input('HOLD')
+
   printMsg('Task PolConvert is Done\n\n')
   return CGains   # RETURN!
 
@@ -3016,9 +3126,14 @@
 
 if __name__=='__main__':
 
-  polconvert(IDI, OUTPUTIDI, DiFXinput, DiFXcalc, doIF, linAntIdx, Range, ALMAant, spw, calAPP, calAPPTime, APPrefant, 
-             gains, interpolation, gainmode,   XYavgTime, dterms, amp_norm, XYadd, XYdel, XYratio, swapXY, swapRL, 
-             IDI_conjugated, plotIF, plotRange, plotAnt,excludeAnts,doSolve,solint,doTest,npix,solveAmp,solveMethod, calstokes, calfield)
-  
-  
-  
+  print('calling polconvert')
+  polconvert(IDI, OUTPUTIDI, DiFXinput, DiFXcalc, doIF, linAntIdx,
+    Range, ALMAant, spw, calAPP, calAPPTime, APPrefant, gains,
+    interpolation, gainmode, XYavgTime, dterms, amp_norm, XYadd, XYdel,
+    XYratio, usePcal, swapXY, swapRL, feedRotation, IDI_conjugated, plotIF,
+    plotRange, plotAnt,excludeAnts,excludeBaselines,doSolve,solint,doTest,
+    npix,solveAmp,solveMethod, calstokes, calfield)
+
+#
+#eof
+#
Index: DataIOSWIN.cpp
===================================================================
--- DataIOSWIN.cpp	(revision 9901)
+++ DataIOSWIN.cpp	(working copy)
@@ -530,8 +530,9 @@
   
   
 // Derive the parallactic angles:
-  getParAng(sidx,ant1-1,ant2-1,UVW,AuxPA1,AuxPA2);
+//z getParAng(sidx,ant1-1,ant2-1,UVW,AuxPA1,AuxPA2);
   daytemp2 = (daytemp + day0)*86400.;
+  getParAng(sidx,ant1-1,ant2-1,UVW,daytemp,AuxPA1,AuxPA2);
 
   
   
Index: DataIO.h
===================================================================
--- DataIO.h	(revision 9901)
+++ DataIO.h	(working copy)
@@ -48,6 +48,8 @@
 double *BaseLine[3];
 double *SinDec;
 double *CosDec;
+//z
+double *RA;
 double *AntLon;
 int *Mount;
 double *Lat;
@@ -94,7 +96,8 @@
    void getFrequencies(double* Freqarray);
 
 // Compute parallactic angle.
-  void getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &P1, double &P2);
+//z  void getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &P1, double &P2);
+   void getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &MJD, double &P1, double &P2);
 
 // Estimate amplitude ratios from the autocorrelations:
   std::complex<float> getAmpRatio(int ant, int spw, int chan);
Index: DataIOFITS.cpp
===================================================================
--- DataIOFITS.cpp	(revision 9901)
+++ DataIOFITS.cpp	(working copy)
@@ -266,7 +266,8 @@
 
  //     printf("READ! \n");
 
-      getParAng(souidx-1,a1-1,a2-1,UVW,AuxPA1,AuxPA2);
+//z    getParAng(souidx-1,a1-1,a2-1,UVW,AuxPA1,AuxPA2);
+      getParAng(souidx-1,a1-1,a2-1,UVW,Times[il],AuxPA1,AuxPA2);
 
  //     printf("COMPUTED: %.2e  %.2e \n",AuxPA1,AuxPA2);
 
@@ -657,7 +658,8 @@
        UVW[0] = (double) FUVW[0]; UVW[1] = (double) FUVW[1]; UVW[2] = (double) FUVW[2]; 
 
 /////// TODO: SORT OUT a1-1 -> a1
-       getParAng(souidx-1,a1-1,a2-1,UVW,AuxPA1,AuxPA2);
+//z    getParAng(souidx-1,a1-1,a2-1,UVW,AuxPA1,AuxPA2);
+       getParAng(souidx-1,a1-1,a2-1,UVW,Times[il],AuxPA1,AuxPA2);
 
     } else if(saveSource<0 || souidx == saveSource) {
 
Index: svn_94.diff
===================================================================
--- svn_94.diff	(nonexistent)
+++ svn_94.diff	(working copy)
@@ -0,0 +1,1010 @@
+Index: _PolConvert.cpp
+===================================================================
+--- _PolConvert.cpp	(revision 9901)
++++ _PolConvert.cpp	(working copy)
+@@ -1,5 +1,4 @@
+ /*
+-
+ # Copyright (c) Ivan Marti-Vidal 2015-2020. 
+ #               EU ALMA Regional Center. Nordic node.
+ #               Observatori Astronomic, Universitat de Valencia
+@@ -52,7 +51,12 @@
+ 
+ */
+  
++//z old code 1, new code 0
++#define IVAN_OLD    1
+ 
++//z straight up merge
++#define IVAN_MERGE  0
++
+ #include <Python.h>
+ 
+ // compiler warning that we use a deprecated NumPy API
+@@ -316,7 +320,8 @@
+   double *AntCoordArr = (double *)PyArray_DATA(antcoordObj);
+   int *AntMountArr = (int *)PyArray_DATA(antmountObj);
+   double *SouCoordArr = (double *)PyArray_DATA(PyList_GetItem(soucoordObj,1));
+-    
++//z
++  double *SouCoordRA = (double *)PyArray_DATA(PyList_GetItem(soucoordObj,0));
+   Geometry->NtotSou = (int) PyArray_DIM(PyList_GetItem(soucoordObj,1),0);
+   Geometry->NtotAnt = (int) PyArray_DIM(antcoordObj,0);
+ 
+@@ -328,6 +333,8 @@
+   Geometry->BaseLine[2] = new double[Nbas+Geometry->NtotAnt+1];
+   Geometry->SinDec = new double[Geometry->NtotSou];
+   Geometry->CosDec = new double[Geometry->NtotSou];
++//z:
++  Geometry->RA = new double[Geometry->NtotSou];
+   Geometry->AntLon = new double[Geometry->NtotAnt];
+   Geometry->Mount = new int[Geometry->NtotAnt];
+   Geometry->Lat = new double[Geometry->NtotAnt];
+@@ -364,8 +371,11 @@
+   };
+ 
+   for (i=0; i<Geometry->NtotSou;i++){
+-      Geometry->SinDec[i] = sin(SouCoordArr[i+Geometry->NtotSou]);
+-      Geometry->CosDec[i] = cos(SouCoordArr[i+Geometry->NtotSou]);
++//z      Geometry->SinDec[i] = sin(SouCoordArr[i+Geometry->NtotSou]);
++//z      Geometry->CosDec[i] = cos(SouCoordArr[i+Geometry->NtotSou]);
++        Geometry->SinDec[i] = sin(SouCoordArr[i]);
++        Geometry->CosDec[i] = cos(SouCoordArr[i]);
++        Geometry->RA[i] = SouCoordRA[i];
+ //      std::cout<< i << " "<< SouCoordArr[i]*180./3.1415926535 << "\n";
+   };
+ 
+@@ -587,13 +597,32 @@
+ 
+ // Which IFs do we convert?
+   int IFs2Conv[nIFconv];
++//z
++#if IVAN_OLD
++#warning "IVAN_OLD CODE selected A"
+   for (ii=0; ii<nIFconv; ii++) {
+-      if (doAll){IFs2Conv[ii]=ii+1;} else {
+-        IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii));
++    if (doAll){IFs2Conv[ii]=ii+1;} else {
++      IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii)); };
++  }
++#else // IVAN_OLD
++#warning "IVAN_NEW CODE selected A"
++  // looks like a shift by 1 here
++  for (ii=0; ii<nIFconv; ii++) {
++    if (doAll){IFs2Conv[ii]=ii;} else {
++      IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii)) - 1; };
++  }
++#endif // IVAN_OLD
++#if IVAN_MERGE
++#warning "IVAN_MERGED CODE selected A"
++  for (ii=0; ii<nIFconv; ii++) {
++//z      if (doAll){IFs2Conv[ii]=ii+1;} else {
++//z        IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii));
++      if (doAll){IFs2Conv[ii]=ii;} else {
++        IFs2Conv[ii] = (int)PyInt_AsLong(PyList_GetItem(doIF,ii)) - 1;
+       };
+   }; 
++#endif // IVAN_MERGE
+ 
+-
+   int nnu = DifXData->getNfreqs();
+   int ALMARefAnt = -1;
+   // If no calAPP is used, do not look for any extra X-Y phase offset.
+@@ -719,6 +748,9 @@
+ 
+ // Prepare plotting files:
+   int noI = -1;
++//z
++#if IVAN_OLD
++#warning "IVAN_OLD CODE selected B"
+   for (ii=0; ii<nIFplot; ii++){
+       sprintf(message,"POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i",IFs2Plot[ii]+1);
+       printf("Writing %s\n", message);
+@@ -729,6 +761,36 @@
+          fwrite(&noI,sizeof(int),1,plotFile[ii]);
+       };
+   };
++#else // IVAN_OLD
++#warning "IVAN_NEW CODE selected B"
++  for (ii=0; ii<nIFconv; ii++) {
++      sprintf(message,"POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i",IFs2Conv[ii]+1);
++      printf("Writing %s\n", message);
++      plotFile[ii] = fopen(message,"wb");
++      if (IFs2Conv[ii]>=0 && IFs2Conv[ii]<nnu){
++         fwrite(&nchans[IFs2Conv[ii]],sizeof(int),1,plotFile[ii]);
++      } else {
++         fwrite(&noI,sizeof(int),1,plotFile[ii]);
++      };
++  };
++#endif // IVAN_OLD
++#if IVAN_MERGE
++#warning "IVAN_MERGED CODE selected B"
++//z  for (ii=0; ii<nIFplot; ii++){
++//z      sprintf(message,"POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i",IFs2Plot[ii]+1);
++  for (ii=0; ii<nIFconv; ii++) {
++      sprintf(message,"POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i",IFs2Conv[ii]+1);
++      printf("Writing %s\n", message);
++      plotFile[ii] = fopen(message,"wb");
++//z      if (IFs2Plot[ii]>=0 && IFs2Plot[ii]<nnu){
++//z         fwrite(&nchans[IFs2Plot[ii]],sizeof(int),1,plotFile[ii]);
++      if (IFs2Conv[ii]>=0 && IFs2Conv[ii]<nnu){
++         fwrite(&nchans[IFs2Conv[ii]],sizeof(int),1,plotFile[ii]);
++      } else {
++         fwrite(&noI,sizeof(int),1,plotFile[ii]);
++      };
++  };
++#endif // IVAN_MERGE
+  
+ 
+   if(doNorm){
+@@ -758,7 +820,20 @@
+ 
+   for (currAntIdx=0; currAntIdx<nALMA; currAntIdx++) {
+       for (im=0; im<nIFconv; im++) {
++//z
++#if IVAN_OLD
++#warning "IVAN_OLD CODE selected C"
+         ii = IFs2Conv[im] - 1;
++#else // IVAN_OLD
++#warning "IVAN_NEW CODE selected C"
++        ii = IFs2Conv[im];
++#endif // IVAN_OLD
++#if IVAN_MERGE
++#warning "IVAN_MERGED CODE selected C"
++//z     ii = IFs2Conv[im] - 1;
++        ii = IFs2Conv[im];
++#endif // IVAN_MERGE
++//z
+         for (ij=0; ij<nchans[ii]; ij++){
+           PrioriGains[currAntIdx][im][ij] *=
+             DifXData->getAmpRatio(currAntIdx, im, ij);
+@@ -780,13 +855,49 @@
+ //  for (ii=0; ii<nnu; ii++) open-brace
+   for (im=0; im<nIFconv; im++) {
+ 
++//z
++#if IVAN_OLD
++#warning "IVAN_OLD CODE selected D"
+     ii = IFs2Conv[im] - 1;
++#else // IVAN_OLD
++#warning "IVAN_NEW CODE selected D"
++    ii = IFs2Conv[im];
++#endif // IVAN_OLD
++#if IVAN_MERGE
++#warning "IVAN_MERGED CODE selected D"
++//z ii = IFs2Conv[im] - 1;
++    ii = IFs2Conv[im];
++#endif // IVAN_MERGE
++//z
+ 
++//z
++#if IVAN_OLD
++#warning "IVAN_OLD CODE selected E"
+     bool plotIF = false;
+     int IFplot = 0;
+     for (ij=0; ij<nIFplot; ij++){
+-      if (IFs2Plot[ij]==ii){plotIF=true;IFplot=ij; break;};
++      if (IFs2Plot[ij]==ii){
++        plotIF=true;
++        IFplot=ij; break;};
+     };
++#else // IVAN_OLD
++#warning "IVAN_NEW CODE selected E"
++    int IFplot = 0;
++    for (ij=0; ij<nIFplot; ij++){
++      if (IFs2Plot[ij]==ii){
++        IFplot=ij; break;};
++#endif // IVAN_OLD
++#if IVAN_MERGE
++#warning "IVAN_MERGED CODE selected E"
++//zz    bool plotIF = false;
++    int IFplot = 0;
++    for (ij=0; ij<nIFplot; ij++){
++      if (IFs2Plot[ij]==ii){
++//zz    plotIF=true;
++        IFplot=ij; break;};
++    };
++#endif // IVAN_MERGE
++//zz
+ 
+  //   if (ii >= nnu) {
+  //     sprintf(message,"ERROR! DATA DO NOT HAVE IF #%i !!\n",ii);  
+@@ -1246,10 +1357,32 @@
+ 
+ // Shall we write in plot file?  auxB -> auxB2 to avoid confusion
+         //indent level within time range
++
++//z
++//z QUESTION: so why does plotIF disappear here?
++//z note auxB -> auxB2
++#if IVAN_OLD
++#warning "IVAN_OLD CODE selected F"
+         auxB2 = (currT>=plRange[0] && currT<=plRange[1] &&
+                 plotIF && (calField<0 || currF==calField));
+                 //plAnt == otherAnt);
++#else // IVAN_OLD
++#warning "IVAN_NEW CODE selected F"
++        auxB2 = (currT>=plRange[0] && currT<=plRange[1] &&
++                (calField<0 || currF==calField));
++                //plAnt == otherAnt);
++#endif // IVAN_OLD
++#if IVAN_MERGE
++#warning "IVAN_MERGED CODE selected F"
++//z     auxB2 = (currT>=plRange[0] && currT<=plRange[1] &&
++//z             plotIF && (calField<0 || currF==calField));
++                //plAnt == otherAnt);
++        auxB2 = (currT>=plRange[0] && currT<=plRange[1] &&
++                (calField<0 || currF==calField));
++                //plAnt == otherAnt);
++#endif // IVAN_MERGE
+ 
++
+         //indent level within time range
+ // Convert:
+         if(Phased){
+Index: _getAntInfo.cpp
+===================================================================
+--- _getAntInfo.cpp	(revision 9901)
++++ _getAntInfo.cpp	(working copy)
+@@ -72,17 +72,25 @@
+     "Returns antenna coordinates";
+ static char getMounts_docstring[] =
+     "Returns the mount types";
++//z
++//static char getNames_docstring[] = 
++//    "Returns the antenna names (codes)";
+ 
+ /* Available functions */
+ static PyObject *getAntInfo(PyObject *self, PyObject *args);
+ static PyObject *getCoords(PyObject *self, PyObject *args);
+ static PyObject *getMounts(PyObject *self, PyObject *args);
++//z
++//static PyObject *getNames(PyObject *self, PyObject *args);
+ 
++
+ /* Module specification */
+ static PyMethodDef module_methods[] = {
+     {"getAntInfo", getAntInfo, METH_VARARGS, getAntInfo_docstring},
+     {"getCoords", getCoords, METH_VARARGS, getCoords_docstring},
+-    {"getMounts", getMounts, METH_VARARGS, getMounts_docstring},
++//z {"getMounts", getMounts, METH_VARARGS, getMounts_docstring},
++    {"getMounts", getMounts, METH_VARARGS, getMounts_docstring}, //,
++//  {"getNames", getNames, METH_VARARGS, getNames_docstring},
+     {NULL, NULL, 0, NULL}   /* terminated by list of NULLs, apparently */
+ };
+ 
+@@ -216,7 +224,8 @@
+   fits_close_file(ifile, &status);
+ 
+   if(status){
+-    printf("\n\nPROBLEM LOSING FITS-IDI!  ERR: %i\n\n",status);
++//z printf("\n\nPROBLEM LOSING FITS-IDI!  ERR: %i\n\n",status);
++    printf("\n\nPROBLEM CLOSING FITS-IDI!  ERR: %i\n\n",status);
+     return Py_BuildValue("i",3);
+   };
+ 
+@@ -226,33 +235,60 @@
+ };
+ 
+ 
++
++
++
++
++
+ static PyObject *getCoords(PyObject *self, PyObject *args){
+ 
+ // Build numpy array:
+ 
+ PyObject *CoordArr;
+-long CD[2] = {Nants,3};
++//z long CD[2] = {Nants,3};
++//z
++//z GBC had already removed Py_INCREF.
++//z Py_INCREF(CoordArr); -- original code
++//z CoordArr = PyArray_SimpleNewFromData(2, &CD[0], NPY_DOUBLE, (void*) Coords);
++//z Py_INCREF(CoordArr); -- needed only if we retain it and use it
++int nd = 2;
++npy_intp* CD = new npy_intp[nd];
+ 
+-//Py_INCREF(CoordArr); -- original code
+-CoordArr = PyArray_SimpleNewFromData(2, &CD[0], NPY_DOUBLE, (void*) Coords);
+-//Py_INCREF(CoordArr); -- needed only if we retain it and use it
++CD[0] = Nants; CD[1] = 3; // = {Nants,3};
+ 
++//Py_INCREF(CoordArr);
++CoordArr = PyArray_SimpleNewFromData(nd, CD, NPY_DOUBLE, (void*) Coords);
++
+ return CoordArr;
+ 
+ };
+ 
+ 
++
++
++
++
++
+ static PyObject *getMounts(PyObject *self, PyObject *args){
+ 
++
+ // Build numpy array:
+ PyObject *MountArr;
+ 
+-long MD[1] = {Nants};
++//z long MD[1] = {Nants};
++int nd = 1;
++npy_intp* MD = new npy_intp[nd];
++MD[0] = Nants;
+ 
+-//Py_INCREF(MountArr); -- original code
+-MountArr = PyArray_SimpleNewFromData(1, &MD[0], NPY_INT, (void*) Mounts);
+-//Py_INCREF(MountArr); -- needed only if we retain it and use it
++//z GBC had already removed Py_INCREF
++//z Py_INCREF(MountArr); -- original code
++//z MountArr = PyArray_SimpleNewFromData(1, &MD[0], NPY_INT, (void*) Mounts);
++//z Py_INCREF(MountArr); -- needed only if we retain it and use it
+ 
++//z
++//Py_INCREF(MountArr);
++MountArr = PyArray_SimpleNewFromData(nd, MD, NPY_INT, (void*) Mounts);
++
+ return MountArr;
+ 
+ };
+Index: DataIO.cpp
+===================================================================
+--- DataIO.cpp	(revision 9901)
++++ DataIO.cpp	(working copy)
+@@ -106,42 +106,63 @@
+ 
+ 
+ 
+-void DataIO::getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &P1, double &P2){
++//z void DataIO::getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &P1, double &P2){
++void DataIO::getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &MJD, double &P1, double &P2){
+ 
+-double V2, Bx, By; //, Bz;
+-double CH, SH, CT1, CT2, HAng, H1, H2, Elev1, Elev2;
+-int ibas;
++//z double V2, Bx, By; //, Bz;
++//z double CH, SH, CT1, CT2, HAng, H1, H2, Elev1, Elev2;
++//zz double V2, Bx, By, GMST; //, Bz;
++double GMST;
++//double CDec, SDec, SRA, TLat1,TLat2, Lon1,Lon2,CH, SH;
++double CT1, CT2, HAng, H1, H2, Elev1, Elev2;
++//zz int ibas;
+ 
+ Elev1 = 0.0; Elev2 = 0.0;
+ 
++//z:
++double days = MJD/86400.;
++double t = (days-51544.0)/36525.;
++double Hh = days - floor(days);
++double GMsec = 24110.54841 + 8640184.812866*t + 0.093104*t*t - 0.0000062*t*t*t;
++GMST = (GMsec/86400. + Hh)*2.*3.1415926535;
++//z:
++//CDec = Geometry->SinDec[sidx];
++//SDec = Geometry->CosDec[sidx];
++//TLat1 = tan(Geometry->Lat[Ant1]);
++//TLat2 = tan(Geometry->Lat[Ant1]);
++//Lon1 = Geometry->AntLon[Ant1];
++//Lon2 = Geometry->AntLon[Ant2];
++//SRA = Geometry->RA[sidx];
+ 
+ // Dummy value if autocorrelation:
+ if (Ant1==Ant2){P1 = -1.e9; P2 = -1.e9; return;};
+ 
++//z more blank lines
+ 
+ if(sidx<Geometry->NtotSou && Ant1<Geometry->NtotAnt && Ant2<Geometry->NtotAnt){
+ 
+-  V2 = Geometry->SinDec[sidx]*UVW[1] - Geometry->CosDec[sidx]*UVW[2];
+-  ibas = Geometry->BasNum[Ant1][Ant2];
++//z  V2 = Geometry->SinDec[sidx]*UVW[1] - Geometry->CosDec[sidx]*UVW[2];
++//zz   ibas = Geometry->BasNum[Ant1][Ant2];
+ 
++//z  if (ibas<0){
++//z  ibas = -ibas;
++//z   Bx = -Geometry->BaseLine[0][ibas];
++//z   By = -Geometry->BaseLine[1][ibas];
++//z// Bz = -Geometry->BaseLine[2][ibas];
++//z  } else {
++//z   Bx = Geometry->BaseLine[0][ibas];
++//z   By = Geometry->BaseLine[1][ibas];
++//z// Bz = Geometry->BaseLine[2][ibas];
++//z  };
++//z
++//z  CH = (UVW[0]*By - V2*Bx); // /(By**2. + Bx**2.);
++//z  SH = (UVW[0]*Bx + V2*By); // /(By**2. + Bx**2.);
+ 
+-  if (ibas<0){
+-  ibas = -ibas;
+-   Bx = -Geometry->BaseLine[0][ibas];
+-   By = -Geometry->BaseLine[1][ibas];
+-//   Bz = -Geometry->BaseLine[2][ibas];
+-  } else {
+-   Bx = Geometry->BaseLine[0][ibas];
+-   By = Geometry->BaseLine[1][ibas];
+-//   Bz = Geometry->BaseLine[2][ibas];
+-  };
+-
+-  CH = (UVW[0]*By - V2*Bx); // /(By**2. + Bx**2.);
+-  SH = (UVW[0]*Bx + V2*By); // /(By**2. + Bx**2.);
+   CT1 = Geometry->CosDec[sidx]*tan(Geometry->Lat[Ant1]);
+   CT2 = Geometry->CosDec[sidx]*tan(Geometry->Lat[Ant2]);
+ 
+-  HAng = atan2(SH,CH);
++//z  HAng = atan2(SH,CH);
++  HAng = GMST - Geometry->RA[sidx];
+   H1 = HAng + Geometry->AntLon[Ant1];
+   H2 = HAng + Geometry->AntLon[Ant2];
+ 
+@@ -148,10 +169,11 @@
+   // sin(lat)*sin(dec)+cos(lat)cos(dec)cos(H)
+ 
+   if (Geometry->Mount[Ant1] > 3){
+-  Elev1 = asin(sin(Geometry->Lat[Ant1])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant1])*Geometry->CosDec[sidx]*cos(H1));};
++   Elev1 = asin(sin(Geometry->Lat[Ant1])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant1])*Geometry->CosDec[sidx]*cos(H1));};
+ 
+   if (Geometry->Mount[Ant2] > 3){
+-  Elev2 = asin(sin(Geometry->Lat[Ant1])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant1])*Geometry->CosDec[sidx]*cos(H1));};
++//z  Elev2 = asin(sin(Geometry->Lat[Ant1])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant1])*Geometry->CosDec[sidx]*cos(H1));};
++   Elev2 = asin(sin(Geometry->Lat[Ant2])*Geometry->SinDec[sidx]+cos(Geometry->Lat[Ant2])*Geometry->CosDec[sidx]*cos(H2));};
+  
+ 
+   switch (Geometry->Mount[Ant1]){
+@@ -178,6 +200,8 @@
+ 
+ //  P2 = atan2(sin(H2), CT2 - Geometry->SinDec[sidx]*cos(H2));
+ 
++//z:
++//if(Ant1==1){printf("%i  %i  %.16e  %.3f  %.3f  %.3f  %.3f\n",Ant1, Ant2, MJD, GMST, HAng, P1, P2);};
+ 
+ } else {
+ 
+Index: _PolGainSolve.cpp
+===================================================================
+--- _PolGainSolve.cpp	(revision 9901)
++++ _PolGainSolve.cpp	(working copy)
+@@ -1654,6 +1654,7 @@
+ 
+        if (af1>=0){
+          RateResVec[af1] += Weights[i][BNum]*(BLRates00[i][BNum] + BLRates11[i][BNum])/2.;
++//z      // printf("%i | %i-%i %.2e  %.2e \n",i, CalAnts[a1],CalAnts[a2], BLRates00[i][BNum],BLRates11[i][BNum]);
+          DelResVec00[af1] += Weights[i][BNum]*BLDelays00[i][BNum];
+          DelResVec11[af1] += Weights[i][BNum]*BLDelays11[i][BNum];
+          Hessian[af1*NantFit + af1] += Weights[i][BNum];
+@@ -1701,21 +1702,36 @@
+ 
+ 
+ printf("\n\nHessian Globalization Matrix:\n\n");
+-bool isSingular, tempSing;
++//zz bool isSingular, tempSing;
++bool isSingular;
+ isSingular=false;
+ for (i=0; i<NantFit; i++){
+   printf("  ");
+-  tempSing = true;
++//z  tempSing = true;
++  if (Hessian[i*NantFit+i]==0.0){isSingular=true;};
+   for (j=0; j<NantFit; j++){
+-     if (Hessian[i*NantFit+j]!=0.0){tempSing=false;};
++//z     if (Hessian[i*NantFit+j]!=0.0){tempSing=false;};
+      printf("%+.2e ",Hessian[i*NantFit+j]);
+   };
+-  if (tempSing){isSingular=true;};
++//z  if (tempSing){isSingular=true;};
++//z printf("\n");
++//z };
+ printf("\n");
+ };
++if(isSingular){printf("Possible singular matrix!\n");};
++printf("\n");
++
++printf("\n\n Residual baseline phase quantities:\n\n");
++for (i=0; i<NantFit; i++){
++  printf("  ");
++  printf(" %.2e | %.2e | %.2e\n",DelResVec00[i],DelResVec11[i],RateResVec[i]);
++};
++
++// GBC added these which may now be redundant
+ if (!isSingular) printf("\n(not Singular)\n");
+ else             printf("\n(is  Singular)\n");
+ 
++printf("\n");
+ 
+ // The Hessian's inverse can be reused for rates, delays and phases!
+ gsl_death_by = GSL_SUCCESS;
+@@ -1753,6 +1769,8 @@
+ }
+ 
+ if(!isSingular){
++//z     printf("Globalizing solutions\n");
++        printf("Globalizing solutions: ");
+         printf("doing gsl_linalg_LU_solve 3x...\n");
+ 	gsl_linalg_LU_solve (&mm.matrix, permm, &RateInd.vector, xx);
+ 	gsl_linalg_LU_solve (&mm.matrix, permm, &Del00Ind.vector, dd0);
+@@ -2274,6 +2292,8 @@
+ 
+   for(k=0;k<NBas;k++){Tm[k]=Times[currIF][0];};
+ 
++//z
++// printf("NVis: %i\n",NVis[currIF]); fflush(stdout);
+ 
+   for (k=0; k<NVis[currIF]; k++){
+ 
+@@ -2436,6 +2456,9 @@
+ 
+     if(BNum>=0){
+ 
++//z
++    //  printf("Fitting for %i %i\n",a1,a2); fflush(stdout);
++
+     AvVis[BNum] += 1;
+ 
+     FeedFactor1 = std::polar(1.0, feedAngle[a1-1])*PA1[currIF][k]; 
+Index: task_polconvert.py
+===================================================================
+--- task_polconvert.py	(revision 9901)
++++ task_polconvert.py	(working copy)
+@@ -51,7 +51,7 @@
+ from __future__ import absolute_import
+ from __future__ import print_function
+ __version__ = "1.8.1"
+-date = 'Oct 11 2020'     
++date = 'Feb 09 2021'     
+ 
+ 
+ ################
+@@ -112,15 +112,24 @@
+ import pickle as pk
+ 
+ if sys.version_info.major < 3:
++  try:
+     from taskinit import *
+     ms = gentools(['ms'])[0]
+     tb = gentools(['tb'])[0]
++  except Exception as ex:
++    print('unable to load casa tools in python2\n\n')
++    raise ex
+ else:
++  try:
+     # from taskutil import *
+     from casatools import ms as ms_casa
+     from casatools import table as tb_casa
+     ms = ms_casa()
+     tb = tb_casa()
++  except Exception as ex:
++    print('unable to load casa tools in python3\n\n')
++    raise ex
++    
+ 
+ #########################################################
+ ###########################
+@@ -134,25 +143,41 @@
+ 
+ if __name__=='__main__':
+ 
++##z not changed yet
++  print('setting up variables')
++
++  # grep name= polconvert.xml | cut -d\" -f 4 | tr \\012 ,
++  global IDI,OUTPUTIDI,DiFXinput,DiFXcalc,doIF,linAntIdx,Range,ALMAant
++  global spw,calAPP,calAPPTime,APPrefant,gains,interpolation,gainmode
++  global XYavgTime,dterms,amp_norm,XYadd,XYdel,XYratio,usePcal,swapXY
++  global swapRL,feedRotation,IDI_conjugated,plotIF,plotRange,plotAnt
++  global excludeAnts,excludeBaselines,doSolve,solint,doTest,npix,solveAmp
++  global solveMethod,calstokes,calfield
++
+   # set things that the task machinery should set 
+   taskname           = "polconvert"
+-  IDI                =  "bm494e-0-b1_1200.difx"
+-  OUTPUTIDI          =  "1200_POLCONVERT.difx"
+-  DiFXinput          =  "bm494e-0-b1_1200.input"
+-  DiFXcalc           =  "bm494e-0-b1_1200.calc"
+-  doIF               =  [35, 36, 37, 38]
++  IDI                =  "no-such-dir.difx"
++  OUTPUTIDI          =  "no-such-dir.difx"
++  DiFXinput          =  "no-such-file.input"
++  DiFXcalc           =  "no-such-file.calc"
++  doIF               =  [6660,6661]
+   linAntIdx          =  [1]
+   Range              =  []
+-  ALMAant            =  "APP_DERIVERABLES/SPECLINE_0.concatenated.ms.ANTENNA"
++  ALMAant            =  "no-such-dir/label.concatenated.ms.ANTENNA"
+   spw                =  0
+-  calAPP             =  "APP_DERIVERABLES/SPECLINE_0.concatenated.ms.calappphase"
++  calAPP             =  "no-such-dir/label.concatenated.ms.calappphase"
+   calAPPTime         =  [0.0, 5.0]
+   APPrefant          =  ""
+-  gains              =  [['APP_DERIVERABLES/SPECLINE_0.concatenated.ms.bandpass-zphs', 'APP_DERIVERABLES/SPECLINE_0.concatenated.ms.flux_inf.APP', 'APP_DERIVERABLES/SPECLINE_0.concatenated.ms.phase_int.APP.XYsmooth', 'APP_DERIVERABLES/SPECLINE_0.calibrated.ms.XY0.APP', 'APP_DERIVERABLES/SPECLINE_0.calibrated.ms.Gxyamp.ALMA']]
++  gains              =  [[
++    'no-such-dir/label.concatenated.ms.bandpass-zphs',
++    'no-such-dir/label.concatenated.ms.flux_inf.APP',
++    'no-such-dir/label.concatenated.ms.phase_int.APP.XYsmooth',
++    'no-such-dir/label.calibrated.ms.XY0.APP',
++    'no-such-dir/label.calibrated.ms.Gxyamp.ALMA']]
+   interpolation      =  []
+   gainmode           =  [['G', 'T', 'G', 'G', 'G']]
+   XYavgTime          =  0.0
+-  dterms             =  ['APP_DERIVERABLES/SPECLINE_0.calibrated.ms.Df0gen.ALMA']
++  dterms             =  ['no-such-dir/label.calibrated.ms.Df0gen.ALMA']
+   amp_norm           =  0.03
+   XYadd              =  {}
+   XYdel              =  {}
+@@ -166,6 +191,7 @@
+   plotRange          =  [0, 0, 0, 0, 2, 0, 0, 0]
+   plotAnt            =  2
+   excludeAnts        =  []
++  excludeBaselines   =  []
+   doSolve            =  -1
+   solint             =  [1, 1]
+   doTest             =  True
+@@ -174,6 +200,22 @@
+   solveMethod        =  "gradient"
+   calstokes          =  [1.0, 0.0, 0.0, 0.0]
+   calfield           =  -1
++
++  try:
++    # rename polconvert.last to polconvert.test to use it
++    pctest = open('polconvert.test')
++    changes = pctest.readlines()
++    pctest.close()
++    exec('\n'.join(changes))
++    print('updated defaults from polconvert.test')
++  except Exception as ex:
++    print('you cannot run this stand-alone without a polconvert.test\n' +
++          'file that re-assigns the variables.  You can use a previous\n' +
++          'polconvert.last if you copy/rename it to polconvert.test\n')
++    raise ex
++
++##z no changes above, yet
++
+ #
+ #
+ #
+@@ -714,15 +756,21 @@
+     printError("ERROR! OUTPUTIDI should be a string!")
+ 
+   if type(plotIF) is int:
+-    plotIF = [plotIF]
++##z
++#   plotIF = [plotIF]
++    if plotIF >0:
++      plotIF = [plotIF]
++    else:
++      plotIF = []
+   for pli in plotIF:
+     if type(pli) is not int:
+       printError("ERROR! plotIF should be an integer or a list of integers!")
+ 
+-  try:
+-    plotAnt = int(plotAnt)
+-  except:
+-    printError("ERROR! plotAnt should be an integer!")
++##z
++# try:
++#   plotAnt = int(plotAnt)
++# except:
++#   printError("ERROR! plotAnt should be an integer!")
+ 
+ 
+   try:
+@@ -857,10 +905,10 @@
+   except:
+     printError("Bad format for calAPPTime. Should be a list of 2 floats!")
+ 
++##z
++# if plotAnt in linAntIdx:
++#   printMsg("WARNING: Plotting will involve autocorrelations. \nThis has not been fully tested!") 
+ 
+-  if plotAnt in linAntIdx:
+-    printMsg("WARNING: Plotting will involve autocorrelations. \nThis has not been fully tested!") 
+-
+   printMsg("XYadd is %s"%str(XYadd))
+   printMsg("XYdel is %s"%str(XYdel))
+   printMsg("XYratio is %s"%str(XYratio))
+@@ -1044,6 +1092,8 @@
+       if decappUnit=='DEGREES':
+         soucoords[1] *= np.pi/180.
+ 
++##z
++      antcodes = [ff[:2] for ff in ffile['ANTENNA'].data['ANNAME']]
+       ffile.close()
+ 
+ # THESE LINES FAIL IF ORBPARM IS PRESENT IN ARRAY GEOMETRY!
+@@ -1057,7 +1107,7 @@
+       else:
+         antcoords = gA.getCoords()
+         antmounts = gA.getMounts()
+-        antcodes = ['%02i'%i for i in range(1,len(antmounts)+1)]
++##z #   antcodes = ['%02i'%i for i in range(1,len(antmounts)+1)]
+     except:
+       printMsg('WARNING! This FITS-IDI file has missing information!\nPolConvert may not calibrate properly.')
+   else:
+@@ -1068,6 +1118,8 @@
+ ######
+ 
+ 
++
++
+ ######
+ # IF THIS IS A SWIN DATASET, READ THE INPUT FILE INTO 
+ # A METADATA LIST:
+@@ -1140,6 +1192,7 @@
+  
+ 
+   else:
++
+ # READ FREQUENCY INFO TO HELP SELECTING THE SPW AUTOMATICALLY:
+     import pyfits as pf
+     fitsf = pf.open(IDI)
+@@ -1149,13 +1202,30 @@
+     IFchan = nch
+     Nr = fitsf['FREQUENCY'].header['NO_BAND']
+     sgn = {True:1.0,False:-1.0}[bw>0.0]
+-    FrInfo = {'FREQ (MHZ)':[nu0/1.e6], 'BW (MHZ)':[bw*nch/1.e6], 'SIGN':[sgn]}
++##z
++##z FrInfo = {'FREQ (MHZ)':[nu0/1.e6], 'BW (MHZ)':[bw*nch/1.e6], 'SIGN':[sgn]}
++    FrInfo = {'FREQ (MHZ)':[], 'BW (MHZ)':[], 'SIGN':[], 'NUM CHANNELS':[]}
++    if sgn:
++      FrInfo['SIDEBAND'] = ['U' for i in range(Nr)]
++    else:
++      FrInfo['SIDEBAND'] = ['L' for i in range(Nr)]
++
++    metadata = []
++##z
+     for i in range(Nr):
+-      FrInfo['FREQ (MHZ)'] += [(nu0 + bw*nch/1.e6)/1.e6]
++##z   FrInfo['FREQ (MHZ)'] += [(nu0 + bw*nch/1.e6)/1.e6]
++      FrInfo['FREQ (MHZ)'] += [(nu0 + i*bw*nch)/1.e6]
+       FrInfo['BW (MHZ)'] += [bw*nch/1.e6]
+       FrInfo['SIGN'] += [sgn]
+       FrInfo['NUM CHANNELS'] += [int(nch)]
++      freqs = nu0 + np.linspace((sgn-1.)/2.,(sgn+1.)/2.,nch,endpoint=False)*bw
++      metadata.append(freqs)
+ 
++    FrInfo['CHANS TO AVG'] = [1 for i in range(Nr)]
++    FrInfo['OVERSAMPLE FAC.'] = [1 for i in range(Nr)]
++    FrInfo['DECIMATION FAC.']=[1 for i in range(Nr)]
++
++
+     if len(doIF)==0:
+      doIF = list(range(1,1+fitsf['FREQUENCY'].header['NO_BAND']))
+ 
+@@ -1163,6 +1233,7 @@
+ 
+ 
+ 
++
+ # ANTENNAS TO PARTICIPATE IN THE GAIN ESTIMATES:
+   nTotAnt = len(antcoords)
+ 
+@@ -1169,16 +1240,41 @@
+   calAnts = []
+   for exA in antcodes:
+     if exA not in excludeAnts:
++##z   #print exA, excludeAnts, exA in excludeAnts
+       calAnts.append(antcodes.index(exA)+1)
+       ### = [i+1 for i in range(len(antcoords)) if i+1 not in excludeAnts]
+     else:
+       printMsg("Excluding antenna %s"%str(exA))
+ 
++##z following section
++  try:
++    plotAnt = int(plotAnt)
++  except:
++    if plotAnt not in antcodes:
++      printError("Reference antenna %s is not found in metadata!"%str(plotAnt))
++    else:
++      plotAnt = antcodes.index(plotAnt)+1
+ 
++  for i in range(len(linAntIdx)):
++    try:
++      linAntIdx[i] = int(linAntIdx[i])
++    except:
++      if linAntIdx[i] not in antcodes:
++        linAntIdx[i] = antcodes.index(linAntIdx[i])+1
++
++
++  if plotAnt in linAntIdx:
++    printMsg(
++        "WARNING: Plotting will involve autocorrelations. \nThis has not been fully tested!")
++
++
++
++
++
+   FlagBas1 = []
+   FlagBas2 = []
+   for fbi in excludeBaselines:
+-    print(fbi, antcodes)
++##z # print(fbi, antcodes)
+     printMsg("Excluding baseline %s"%str(fbi))
+     if fbi[0] in antcodes and fbi[1] in antcodes:
+       FlagBas1.append(antcodes.index(fbi[0])+1) ### = np.array([int(i[0]+1) for i in excludeBaselines])
+@@ -1211,7 +1307,6 @@
+ 
+ 
+ 
+-
+ #######################
+ ##### GET SPECTRAL WINDOW AUTOMATICALLY:
+   if isPhased and spw < 0:
+@@ -1401,6 +1496,7 @@
+    except:
+      printError("Bad time range format for plotRange!")
+ 
++
+   if len(Range) == 0:
+     Ran = np.array([0.,1.e20])
+   else:
+@@ -1910,7 +2006,7 @@
+ 
+ 
+   if DEBUG:
+-    PC_Params = [nALMATrue, doIF, plotAnt, len(allants), doIF,
++    PC_Params = [nALMATrue, plotIF, plotAnt, len(allants), doIF,
+         swapXY, ngain, NSUM, kind, len(gaindata), len(dtdata), OUTPUT,
+         linAntIdxTrue, plRan, Ran, allantidx, len(nphtimes), len(antimes),
+         len(refants), len(asdmtimes),  doTest, doSolve, doConj, doAmpNorm,
+@@ -1935,6 +2031,7 @@
+ 
+   printMsg("\n###\n### Done with PolConvert (status %d).\n###" % (didit))
+ 
++##z # raw_input('HOLD')
+ 
+   if didit != 0:
+     printError("\n\n ERROR IN POLCONVERT!\n\n")
+@@ -2136,8 +2233,8 @@
+ 
+        i += 1
+ 
+-       while Chi2_0>currChi2: 
+-
++##z    while Chi2_0>currChi2: 
++       while Chi2_0>=currChi2:
+          i += 1     
+          LMTune *= KFacRaise
+          ptst0 = np.copy(currP)
+@@ -2148,8 +2245,13 @@
+          if i>=MAXIT:
+            break
+ 
+-       relchange = (currChi2 - Chi2_0)/Chi2_0
++##z    relchange = (currChi2 - Chi2_0)/Chi2_0
++       if Chi2_0 >0.0:
++         relchange = (currChi2 - Chi2_0)/Chi2_0
++       else:
++         printError("\n\n  Problem in PolGainSolve.\n")         
+ 
++
+ # No improvement:
+        if currChi2<minChi2:
+          minChi2 = currChi2
+@@ -2217,7 +2319,10 @@
+ 
+      Chi2_final = PS.GetChi2(minGains,LMTune,Ch0,Ch1,1)
+ 
+-     FLIP = Chi2_final < 0.0  # Flip gains by 180 degrees.
++##z  FLIP = Chi2_final < 0.0  # Flip gains by 180 degrees.
++     FLIP = False # Chi2_final < 0.0  # Flip gains by 180 degrees.
++
++# GBC debugging:
+      if FLIP: sys.stdout.write(' Flipped\n')
+      else:    sys.stdout.write(' NotFlip\n')
+      printMsg("    Final error: %.3e in ChSq"%(np.abs(relchange)))
+@@ -2456,7 +2561,8 @@
+       sub2.set_xlim((np.min(Freq2Plot) - Dnu*0.1,np.max(Freq2Plot) + Dnu*0.45))
+       sub2.set_ylim((0.,2.5))
+ 
+-      sub2.legend(numpoints=1)
++##z   sub2.legend(numpoints=1)
++      sub1.legend(numpoints=1)
+       sub1.set_ylabel('Cross-Phase (deg.)')
+       sub2.set_ylabel('Cross-Amp (Norm.)')
+       sub2.set_xlabel('Frequency (GHz)')
+@@ -2536,7 +2642,8 @@
+      if os.stat("POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i"%pli).st_size>10:
+        GoodIFs.append(pli)
+      else:
+-       printMsg("WARNING! IF %i was NOT polconverted properly\n"%pli)
++       printMsg(("WARNING! IF %i was NOT polconverted properly\n"%pli) +
++         ("POLCONVERT.FRINGE/POLCONVERT.FRINGE_%i missing"%pli))
+ 
+ # start of GoodIFs pli loop
+    for pli in GoodIFs:
+@@ -3010,6 +3117,9 @@
+     printMsg(str(ex))
+   ofile.close()
+   printMsg('PolConvert.XYGains.dat was written with CGains' + str(CGains.keys()))
++
++# raw_input('HOLD')
++
+   printMsg('Task PolConvert is Done\n\n')
+   return CGains   # RETURN!
+ 
+@@ -3016,9 +3126,14 @@
+ 
+ if __name__=='__main__':
+ 
+-  polconvert(IDI, OUTPUTIDI, DiFXinput, DiFXcalc, doIF, linAntIdx, Range, ALMAant, spw, calAPP, calAPPTime, APPrefant, 
+-             gains, interpolation, gainmode,   XYavgTime, dterms, amp_norm, XYadd, XYdel, XYratio, swapXY, swapRL, 
+-             IDI_conjugated, plotIF, plotRange, plotAnt,excludeAnts,doSolve,solint,doTest,npix,solveAmp,solveMethod, calstokes, calfield)
+-  
+-  
+-  
++  print('calling polconvert')
++  polconvert(IDI, OUTPUTIDI, DiFXinput, DiFXcalc, doIF, linAntIdx,
++    Range, ALMAant, spw, calAPP, calAPPTime, APPrefant, gains,
++    interpolation, gainmode, XYavgTime, dterms, amp_norm, XYadd, XYdel,
++    XYratio, usePcal, swapXY, swapRL, feedRotation, IDI_conjugated, plotIF,
++    plotRange, plotAnt,excludeAnts,excludeBaselines,doSolve,solint,doTest,
++    npix,solveAmp,solveMethod, calstokes, calfield)
++
++#
++#eof
++#
+Index: DataIOSWIN.cpp
+===================================================================
+--- DataIOSWIN.cpp	(revision 9901)
++++ DataIOSWIN.cpp	(working copy)
+@@ -530,8 +530,9 @@
+   
+   
+ // Derive the parallactic angles:
+-  getParAng(sidx,ant1-1,ant2-1,UVW,AuxPA1,AuxPA2);
++//z getParAng(sidx,ant1-1,ant2-1,UVW,AuxPA1,AuxPA2);
+   daytemp2 = (daytemp + day0)*86400.;
++  getParAng(sidx,ant1-1,ant2-1,UVW,daytemp,AuxPA1,AuxPA2);
+ 
+   
+   
+Index: DataIO.h
+===================================================================
+--- DataIO.h	(revision 9901)
++++ DataIO.h	(working copy)
+@@ -48,6 +48,8 @@
+ double *BaseLine[3];
+ double *SinDec;
+ double *CosDec;
++//z
++double *RA;
+ double *AntLon;
+ int *Mount;
+ double *Lat;
+@@ -94,7 +96,8 @@
+    void getFrequencies(double* Freqarray);
+ 
+ // Compute parallactic angle.
+-  void getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &P1, double &P2);
++//z  void getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &P1, double &P2);
++   void getParAng(int sidx, int Ant1, int Ant2, double*UVW, double &MJD, double &P1, double &P2);
+ 
+ // Estimate amplitude ratios from the autocorrelations:
+   std::complex<float> getAmpRatio(int ant, int spw, int chan);
+Index: DataIOFITS.cpp
+===================================================================
+--- DataIOFITS.cpp	(revision 9901)
++++ DataIOFITS.cpp	(working copy)
+@@ -266,7 +266,8 @@
+ 
+  //     printf("READ! \n");
+ 
+-      getParAng(souidx-1,a1-1,a2-1,UVW,AuxPA1,AuxPA2);
++//z    getParAng(souidx-1,a1-1,a2-1,UVW,AuxPA1,AuxPA2);
++      getParAng(souidx-1,a1-1,a2-1,UVW,Times[il],AuxPA1,AuxPA2);
+ 
+  //     printf("COMPUTED: %.2e  %.2e \n",AuxPA1,AuxPA2);
+ 
+@@ -657,7 +658,8 @@
+        UVW[0] = (double) FUVW[0]; UVW[1] = (double) FUVW[1]; UVW[2] = (double) FUVW[2]; 
+ 
+ /////// TODO: SORT OUT a1-1 -> a1
+-       getParAng(souidx-1,a1-1,a2-1,UVW,AuxPA1,AuxPA2);
++//z    getParAng(souidx-1,a1-1,a2-1,UVW,AuxPA1,AuxPA2);
++       getParAng(souidx-1,a1-1,a2-1,UVW,Times[il],AuxPA1,AuxPA2);
+ 
+     } else if(saveSource<0 || souidx == saveSource) {
+