diff --git a/contrib/stack/stripmapStack/Stack.py b/contrib/stack/stripmapStack/Stack.py
index 276494e7b..fe55e661a 100755
--- a/contrib/stack/stripmapStack/Stack.py
+++ b/contrib/stack/stripmapStack/Stack.py
@@ -474,19 +474,21 @@ def denseOffsets_Network(self, pairs, stackReference, secondaryDates, config_pre
             self.runf.write(self.text_cmd + 'stripmapWrapper.py -c '+ configName+'\n')
 
   
-    def invertMisregPoly(self):
+    def invertMisregPoly(self, stackReference):
 
         pairDirs = os.path.join(self.workDir, 'refineSecondaryTiming/pairs/')
         dateDirs = os.path.join(self.workDir, 'refineSecondaryTiming/dates/')
-        cmd = self.text_cmd + 'invertMisreg.py -i ' + pairDirs + ' -o ' + dateDirs
+        cmd = self.text_cmd + 'invertMisreg.py -i ' + pairDirs + ' -o ' + dateDirs + ' -r ' + stackReference
         self.runf.write(cmd + '\n')
 
 
-    def  invertDenseOffsets(self):
+    def  invertDenseOffsets(self, stackReference=None):
 
         pairDirs = os.path.join(self.workDir, self.dense_offsets_folder, 'pairs')
         dateDirs = os.path.join(self.workDir, self.dense_offsets_folder, 'dates')
         cmd = self.text_cmd + 'invertOffsets.py -i ' + pairDirs + ' -o ' + dateDirs
+        if stackReference is not None:
+            cmd += ' -r ' + stackReference
         self.runf.write(cmd + '\n')
 
     def rubbersheet(self, secondaryDates, config_prefix):
diff --git a/contrib/stack/stripmapStack/invertMisreg.py b/contrib/stack/stripmapStack/invertMisreg.py
index 8be7e9336..e73e069c6 100755
--- a/contrib/stack/stripmapStack/invertMisreg.py
+++ b/contrib/stack/stripmapStack/invertMisreg.py
@@ -24,6 +24,8 @@ def createParser():
             help='Directory with the overlap directories that has calculated misregistration for each pair')
     parser.add_argument('-o', '--output', type=str, dest='output', required=True,
             help='output directory to save misregistration for each date with respect to the stack Reference date')
+    parser.add_argument('-r', '--reference', type=str, dest='reference', default=None,
+            help='Reference date (YYYYMMDD format). If not provided, the first date in sorted dateList will be used.')
 #    parser.add_argument('-f', '--misregFileName', type=str, dest='misregFileName', default='misreg.txt',
 #            help='misreg file name that contains the calculated misregistration for a pair')
 
@@ -93,7 +95,7 @@ def getPolyInfo(filename):
   #return np.size(azCoefs), np.size(rgCoefs), np.shape(azCoefs), np.shape(rgCoefs)
 
 ######################################
-def design_matrix(pairDirs):
+def design_matrix(pairDirs, referenceDate=None):
   '''Make the design matrix for the inversion.  '''
   tbase,dateList,dateDict = date_list(pairDirs)
   numDates = len(dateDict)
@@ -125,12 +127,20 @@ def design_matrix(pairDirs):
     Laz[ni,:] = azOff[:]
     Lrg[ni,:] = rgOff[:]
 
-  A = A[:,1:]
-  B = B[:,:-1]
+  if referenceDate is None:
+    refIdx = 0
+  elif referenceDate not in dateList:
+    print('Warning: Reference date {} not found in dateList. Using first date as reference.'.format(referenceDate))
+    refIdx = 0
+  else:
+    refIdx = dateList.index(referenceDate)
+
+  A = np.delete(A, refIdx, axis=1)
+  B = np.delete(B, refIdx, axis=1)
   
  # ind=~np.isnan(Laz)
  # return A[ind[:,0],:],B[ind[:,0],:],Laz[ind,:], Lrg[ind]
-  return A, B, Laz, Lrg
+  return A, B, Laz, Lrg, refIdx
  
 ######################################
 def main(iargs=None):
@@ -143,7 +153,11 @@ def main(iargs=None):
 
   tbase, dateList, dateDict = date_list(pairDirs)
 
-  A, B, Laz, Lrg = design_matrix(pairDirs)
+  referenceDate = inps.reference
+  if referenceDate is None:
+    referenceDate = dateList[0]
+
+  A, B, Laz, Lrg, refIdx = design_matrix(pairDirs, referenceDate=referenceDate)
   A1 = np.linalg.pinv(A)
   A1 = np.array(A1,np.float32)
 
@@ -158,8 +172,10 @@ def main(iargs=None):
   RMSE_az = np.sqrt(np.sum(residual_az**2)/len(residual_az))
   RMSE_rg = np.sqrt(np.sum(residual_rg**2)/len(residual_rg))
 
-  Saz = np.vstack((np.zeros((1,Saz.shape[1]), dtype=np.float32), Saz))
-  Srg = np.vstack((np.zeros((1,Srg.shape[1]), dtype=np.float32), Srg))
+  zero_row_az = np.zeros((1, Saz.shape[1]), dtype=np.float32)
+  zero_row_rg = np.zeros((1, Srg.shape[1]), dtype=np.float32)
+  Saz = np.vstack([Saz[:refIdx], zero_row_az, Saz[refIdx:]])
+  Srg = np.vstack([Srg[:refIdx], zero_row_rg, Srg[refIdx:]])
 
   print('')
   print('Rank of design matrix: ' + str(np.linalg.matrix_rank(A)))
diff --git a/contrib/stack/stripmapStack/invertOffsets.py b/contrib/stack/stripmapStack/invertOffsets.py
index 9f13221a8..2c769cb5d 100755
--- a/contrib/stack/stripmapStack/invertOffsets.py
+++ b/contrib/stack/stripmapStack/invertOffsets.py
@@ -29,6 +29,8 @@ def createParser():
             help='Directory with the pair directories that includes dense offsets for each pair')
     parser.add_argument('-o', '--output', type=str, dest='output', required=True,
             help='output directory to save dense-offsets for each date with respect to the stack Reference date')
+    parser.add_argument('-r', '--reference', type=str, dest='reference', default=None,
+            help='Reference date (YYYY-MM-DD HH:MM:SS format from h5 file, or YYYYMMDD format). If not provided, the first date in sorted dateList will be used.')
 
     return parser
 
@@ -104,7 +106,7 @@ def date_list(h5file):
 
 #####################################
 
-def design_matrix(h5File):
+def design_matrix(h5File, referenceDate=None):
   tbase,dateList,dateDict, references, secondarys = date_list(h5File)
   numDates = len(dateDict)
   numPairs = len(references)
@@ -118,17 +120,31 @@ def design_matrix(h5File):
      A[ni,ndxt2] = 1
      B[ni,ndxt1:ndxt2] = tbase[ndxt1+1:ndxt2+1]-tbase[ndxt1:ndxt2]
 
-  #print('A',A)
-  #print('%%%%%%%%%%%%%%% %%%%%')  
-  A = A[:,1:]
-  B = B[:,:-1]
+  if referenceDate is not None and referenceDate not in dateList:
+    referenceDate = datetime.datetime(
+        *time.strptime(referenceDate, "%Y%m%d")[0:5]
+    ).strftime("%Y-%m-%d %H:%M:%S")
 
-  return A, B  
+  if referenceDate is None:
+    refIdx = 0
+  elif referenceDate not in dateList:
+    print('Warning: Reference date {} not found in dateList. Using first date as reference.'.format(referenceDate))
+    refIdx = 0
+  else:
+    refIdx = dateList.index(referenceDate)
+
+  A = np.delete(A, refIdx, axis=1)
+  B = np.delete(B, refIdx, axis=1)
+
+  return A, B, refIdx  
 
 def invert_wlq(inps,h5File):
     tbase,dateList,dateDict, references, secondarys = date_list(h5File)
+    referenceDate = inps.reference
+    if referenceDate is None:
+      referenceDate = dateList[0]
     numPairs = len(references)
-    A,B = design_matrix(h5File)
+    A,B,refIdx = design_matrix(h5File, referenceDate=referenceDate)
    
     h5 = h5py.File(h5File,'r')
     data = h5['/platform-track/observations'].get('offset-azimuth')
@@ -168,9 +184,12 @@ def invert_wlq(inps,h5File):
         #Cm = np.vstack((np.zeros((1,ts.shape[1]), dtype=np.float32), ts))
 
            ts = ts.reshape([NumValidPixels,Npar]).T
-        #ts = np.vstack((np.zeros((1,ts.shape[1]), dtype=np.float32), ts))
-           ds[1:,j,ind] = ts
-           dsq[1:,j,ind] = Cm
+           zero_row = np.zeros((1, ts.shape[1]), dtype=np.float32)
+           ts = np.vstack([ts[:refIdx], zero_row, ts[refIdx:]])
+           ds[:,j,ind] = ts
+           Cm_zero = np.zeros((1, Cm.shape[1]), dtype=np.float32)
+           Cm = np.vstack([Cm[:refIdx], Cm_zero, Cm[refIdx:]])
+           dsq[:,j,ind] = Cm
 
     dateListE = [d.encode("ascii", "ignore") for d in dateList]
     dateListE = np.array(dateListE)
@@ -219,8 +238,11 @@ def invert_wlq(inps,h5File):
 def invert(inps,h5File):
 
     tbase,dateList,dateDict, references, secondarys = date_list(h5File)
+    referenceDate = inps.reference
+    if referenceDate is None:
+      referenceDate = dateList[0]
     numPairs = len(references)
-    A,B = design_matrix(h5File)
+    A,B,refIdx = design_matrix(h5File, referenceDate=referenceDate)
 
     h5 = h5py.File(h5File,'r')
     data = h5['/platform-track/observations'].get('offset-azimuth')    
@@ -245,7 +267,8 @@ def invert(inps,h5File):
         #dsr[:,i,:] =  L_residual
         dst[i,:] = np.absolute(np.sum(np.exp(1j*L_residual),0))/Nz
 
-        ts = np.vstack((np.zeros((1,ts.shape[1]), dtype=np.float32), ts))
+        zero_row = np.zeros((1, ts.shape[1]), dtype=np.float32)
+        ts = np.vstack([ts[:refIdx], zero_row, ts[refIdx:]])
         ds[:,i,:] = ts
 
     dateListE = [d.encode("ascii", "ignore") for d in dateList]
diff --git a/contrib/stack/stripmapStack/stackStripMap.py b/contrib/stack/stripmapStack/stackStripMap.py
index 6a57dd75e..f06084bde 100755
--- a/contrib/stack/stripmapStack/stackStripMap.py
+++ b/contrib/stack/stripmapStack/stackStripMap.py
@@ -182,7 +182,7 @@ def slcStack(inps, acquisitionDates, stackReferenceDate, secondaryDates, pairs,
     i+=1
     runObj = run()
     runObj.configure(inps, 'run_{:02d}_invertMisreg'.format(i))
-    runObj.invertMisregPoly()
+    runObj.invertMisregPoly(stackReferenceDate)
     runObj.finalize()
 
     i+=1
@@ -203,7 +203,7 @@ def slcStack(inps, acquisitionDates, stackReferenceDate, secondaryDates, pairs,
         i+=1
         runObj = run()
         runObj.configure(inps, 'run_{:02d}_invertDenseOffsets'.format(i))
-        runObj.invertDenseOffsets()
+        runObj.invertDenseOffsets(stackReferenceDate)
         runObj.finalize()
 
         i+=1