Import of existing CoSAR code

.gitignore

+
+*.pyc

cosar/__init__.py

+"""
+===================================
+Input and output (:mod:`pycosar`)
+===================================
+
+.. currentmodule:: pycosar
+
+   
+"""
+
+from simplesim import *
+from results import *
+from geosim1d import *
+#from rat import *
+
+__all__ = ['simplesim','results']
\ No newline at end of file

cosar/geosim1d.py

+# -*- coding: utf-8 -*-
+"""
+Created on Wed Dec 18 16:13:58 2013
+
+@author: lope_fr
+"""
+
+
+import numpy as np
+import matplotlib
+import matplotlib.pyplot as plt
+from scipy import constants
+from drama import utils
+# csl stands fro cosarlib
+import cosar.simplesim as csl
+import matplotlib
+matplotlib.rcParams['ps.fonttype'] = 42
+
+def geosim1d(cosar, surface, dx, n_rg=1, PRF=0, Tint=300.0, hamming=True,
+             n_lags=1, surfaceseed=None, compute_IRF=False, pos_pt=[]):
+    """
+    Simple 1D simulation using the geometry defined in cosar
+    @type cosar csl.CoSARGeometry
+    """
+    #Init cosar if nothing given. The secquence also forces the
+    #assert isinstance(cosar, csl.CoSARGeometry)
+    PRF_min = 2*cosar.dDoppler_bandwidth()
+    if PRF == 0:
+        PRF = PRF_min
+    else:
+        if PRF < PRF_min:
+            print("PRF given to low")
+            PRF = PRF_min
+    #Check that the azimuth resolution is finer than the minimum lenth scale
+    if dx > surface.L_min()/2:
+        dx = surface.L_min()/2
+    n_az = int(np.round(cosar.L_illuminated()/dx))
+    n_t = int(Tint * PRF)
+    msg = "Number of azimuth grid samples: " + repr(n_az)
+    print(msg)
+    msg = "Number of pulses: " + repr(n_t)
+    print(msg)
+    k0 = 2*np.pi*cosar.f0/constants.c
+    #Seed for repeatibility of the surface
+    np.random.seed(surfaceseed)
+    nrcs, dop, h, az = surface.scene(cosar.L_illuminated(), dx, cosar.f0)
+    if compute_IRF:
+        nrcs[:] = 0
+        if np.size(pos_pt) == 0:
+            pos_pt = [0]
+        #nrcs[n_az/2 + np.round(1e3/dx)] = 1
+        dop[:] = 0
+        h[:] = 0
+    for pos_this_pt in pos_pt:
+        pos_this_pt_s = np.round(pos_this_pt/dx) + n_az/2
+        nrcs[pos_this_pt_s] = 1
+    #Reseed to randomize realization
+    np.random.seed()
+    #Time vector
+    t_v = np.linspace(-Tint / 2, Tint / 2, int(Tint * PRF))
+    #Doppler phase shift for each point on the ground
+    dop_phasor = np.exp(2j * np.pi * dop / PRF)
+    #The expected amplitude takes into account the grid size and the fact
+    #that we are working with circular Gaussian signals
+    expected_amp = np.sqrt(nrcs * dx / 2)
+    #Now we calculate the trajectory of the two spacecraft assuming linear
+    #motions
+    pos1 = np.zeros([n_t, 3])
+    pos2 = np.zeros([n_t, 3])
+    pos1[:, 1] = -np.sin(cosar.theta_i) * cosar.R0()
+    pos2[:, 1] = -np.sin(cosar.theta_i - cosar.dtheta_i()) * cosar.R0()
+    pos1[:, 2] = np.cos(cosar.theta_i) * cosar.R0()
+    pos2[:, 2] = np.cos(cosar.theta_i - cosar.dtheta_i()) * cosar.R0()
+    for dim_ind in range(3):
+        pos1[:, dim_ind] = pos1[:, dim_ind] + cosar.dv[dim_ind] / 2.0 * t_v
+        pos2[:, dim_ind] = pos2[:, dim_ind] - cosar.dv[dim_ind] / 2.0 * t_v
+    #Correlation coefficcient
+    # alpha**(tau*PRF) = 1/e
+    # tau*PRF * log(alpha)=-1 alpha = exp(-1/(tau*PRF)
+    alpha = 1.0
+    if surface.tau > 0:
+        alpha = np.exp(-1.0/(surface.tau*PRF))
+        msg = "Pulse to pulse coherence: " + repr(alpha)
+        print(msg)
+        alpha_end = np.exp(-1.0*t_v.size/(surface.tau*PRF))
+        msg = "Start to end coherence: " + repr(alpha_end)
+        print(msg)
+    if n_rg > 1:
+        uscat = np.random.randn(n_az, n_rg) + 1j*np.random.randn(n_az, n_rg)
+    else:
+        uscat = np.sqrt(2) * np.exp(1j*2*np.pi*np.random.rand(n_az, n_rg))
+            #Now main loop
+    for t_ind in range(0, t_v.size):
+        #Generate random complex scatterers
+        uscat = (alpha * uscat * dop_phasor.reshape(n_az, 1) +
+                 np.sqrt(1-alpha**2) * (np.random.randn(n_az, n_rg) +
+                                        1j * np.random.randn(n_az, n_rg)))
+        #Multipliy them with a fixed amplitude
+        scat = uscat * expected_amp.reshape(n_az, 1)
+        #Calculate ranges from targets to both radars
+        r1 = ((np.sqrt((pos1[t_ind, 1])**2 +
+              (az - pos1[t_ind, 0])**2 +
+              (pos1[t_ind, 2] - h)**2)).reshape(n_az, 1))
+        r2 = ((np.sqrt((pos2[t_ind, 1])**2 +
+              (az - pos2[t_ind, 0])**2 +
+              (pos2[t_ind, 2] - h)**2)).reshape(n_az, 1))
+        #Include 2-way phase in field
+        s1 = scat * np.exp(-2j*k0*r1)
+        s2 = scat * np.exp(-2j*k0*r2)
+        # Compute range profiles by integrating in azimuth.
+        # This assumed that there is no range migration
+        rx1_ = s1.sum(axis=0)
+        rx2_ = s2.sum(axis=0)
+        # Save in a temporal stack
+        if t_ind == 0:
+            rx1 = rx1_.copy()
+            rx2 = rx2_.copy()
+        else:
+            rx1 = np.vstack((rx1, rx1_))
+            rx2 = np.vstack((rx2, rx2_))
+    #Now we have a range-slow time stack of CoSAR data
+    n_az_out = n_az
+    proc_img = np.zeros([2 * n_lags + 1, n_az_out], dtype=np.complex)
+    for lag_ind in range(0, 2*n_lags + 1):
+        lag = lag_ind - n_lags
+        proc_img[lag_ind, :], T_i_scaling = \
+            csl.cosar_proc_1d(rx1, rx2, pos1, pos2, az, k0,
+                              han=hamming, dt=lag)
+    #Backprojection SAR focusing of one channel
+    sar_img, kk = csl.sar_proc_1d(rx1, pos1, az, k0, han=hamming)
+    #Amplitude scaling of processed image
+    C_p = 2 * cosar.dv[0] / (cosar.wavelength() * cosar.R0())
+    amp_scaling = C_p / PRF
+    proc_img = amp_scaling * proc_img
+    sar_img = amp_scaling * sar_img
+    T_i_scaled = T_i_scaling * Tint
+    return (az, nrcs, dop, h, proc_img, n_lags, T_i_scaled, dx, n_rg, sar_img)
+
+
+def geosim1d_montecarlo(cosar, surface, dx, n_rg=1, PRF=0, Tint=300.0,
+                        hamming=True, n_lags=1, n_it=10, surfaceseed=1):
+    """
+    A function to make several runs
+    """
+    sim_res = geosim1d(cosar, surface, dx, n_rg=n_rg, PRF=PRF, Tint=Tint,
+                       n_lags=n_lags, hamming=hamming, surfaceseed=surfaceseed)
+    #Unpack first simulation
+    az, nrcs, dop, h, proc_img, n_lags, T_i, dx, n_rg, sar_img = sim_res
+    dim_res = proc_img.shape
+    dim_mtc = (n_it,) + dim_res
+    proc_imgs = np.zeros(dim_mtc, dtype=np.complex)
+    #Insert output of first simulation in stack
+    proc_imgs[0, :, :] = proc_img
+    for it in range(1, n_it):
+        sim_res = geosim1d(cosar, surface, dx, n_rg=n_rg, PRF=PRF, Tint=Tint,
+                           n_lags=n_lags, hamming=hamming,
+                           surfaceseed=surfaceseed)
+        az, nrcs, dop, h, proc_img, n_lags, T_i, dx, n_rg, sar_img = sim_res
+        proc_imgs[it, :, :] = proc_img
+
+    return (az, nrcs, dop, h, proc_imgs, n_lags, T_i, dx, n_rg, PRF)
+
+
+def geosim1d_save(sim_results, fileroot):
+    """
+    Function to save simulator results
+    """
+    #Unpack simulation results
+    az, nrcs, dop, h, proc_imgs, n_lags, T_i, dx, n_rg, PRF = sim_results
+    filename = fileroot+'_res.npz'
+    np.savez(filename, az=az, nrcs=nrcs, dop=dop, h=h, proc_imgs=proc_imgs,
+             n_lags=n_lags, T_i=T_i, dx=dx, n_rg=n_rg, PRF=PRF)
+
+
+def geosim1d_load(fileroot):
+    """
+    Function to save simulator results
+    """
+    filename = fileroot+'_res.npz'
+    res_d = np.load(filename)
+    res = (res_d['az'], res_d['nrcs'], res_d['dop'], res_d['h'],
+           res_d['proc_imgs'], res_d['n_lags'], res_d['T_i'],
+           res_d['dx'], res_d['n_rg'], res_d['PRF'])
+    return (res)
+
+
+def geosim1d_analyze(cosar, surface, sim_results, dB=False, fileroot=False,
+                     fontsize=16):
+    """
+    Function to analyze output of geosim1d
+    """
+    matplotlib.rcParams.update({'font.size': fontsize})
+    #Unpack simulation results
+    az, nrcs, dop, h, proc_imgs, n_lags, T_i, dx, n_rg, PRF = sim_results
+    #Calculate Quality factor
+    az_res = 0.5 * cosar.wavelength() * cosar.R0() / (cosar.dv[0] * T_i)
+    #tau_ca
+    tau_ca_Doppler = (cosar.wavelength()*np.sqrt(2) /
+                      (4 * np.pi * surface.vr_stdev))
+    tau_ca = np.sqrt(1.0 / (1 / surface.tau**2 + 1 / tau_ca_Doppler**2))
+    tau_ca_quant = np.max([tau_ca, 1 / PRF])
+    Q_nrcs = csl.Q_w(nrcs, nrcs, dx, n_rg, T_i, tau_ca_quant, az_res)
+    #Expected standard deviation of the estimated nrcs
+    nrcs_err_t = nrcs / np.sqrt(Q_nrcs)
+    #Lag 1 expected value
+    dop_phase = 2 * np.pi * dop / PRF
+    dop_phasor = np.exp(2j * np.pi * dop / PRF)
+    #Height phasor
+    h_phase = cosar.kz() * h
+    h_phasor = np.exp(1j * h_phase)
+    corr_ati = 1.0
+    if surface.tau > 0:
+        corr_ati = np.exp(-1.0 / (surface.tau * PRF))
+    R_tau_lag1 = corr_ati * dop_phasor * h_phasor * nrcs
+    R_tau_lag1m = corr_ati * np.conj(dop_phasor) * h_phasor * nrcs
+    R_tau_lag0 = h_phasor * nrcs
+    if np.size(proc_imgs.shape) == 3:
+        proc_img = proc_imgs[0, :, :]
+        nrcs_err = proc_imgs[:, n_lags, :] - R_tau_lag0.reshape((1, nrcs.size))
+        nrcs_err_mean = np.mean(nrcs_err, axis=0)
+        nrcs_err_std = np.std(nrcs_err, axis=0)
+        plt.figure()
+        ax = plt.subplot(1, 1, 1)
+        p_std, = plt.plot(az / 1e3, nrcs_err_std, 'b',
+                          label=r'std($\delta R_{\tau}$)')
+        p_mean_r, = plt.plot(az / 1e3, np.real(nrcs_err_mean), 'c',
+                             label=r'$\mathrm{Re}\{\overline{\delta R_{\tau}}\}$')
+        p_mean_i, = plt.plot(az / 1e3, np.imag(nrcs_err_mean), 'm',
+                             label=r'$\mathrm{Im}\{\overline{\delta R_{\tau}}\}$')
+        plt.plot(az / 1e3, nrcs_err_t, 'b--')
+        plt.xlabel('Azimuth [km]', fontsize=fontsize)
+        plt.ylabel('NRCS errors', fontsize=fontsize)
+        ax.set_ylim([-0.05,0.1])
+        handles, labels = ax.get_legend_handles_labels()
+        l1 = plt.legend([handles[0]], [labels[0]], loc=2)
+        l2 = plt.legend(handles[1:], labels[1:], loc=3, ncol=2)
+        plt.gca().add_artist(l1)
+        if type(fileroot) == str:
+            filename = fileroot+'_err_stats.eps'
+            plt.savefig(filename)
+            filename = fileroot+'_err_stats.png'
+            plt.savefig(filename)
+        plt.figure()
+        plt.plot(az/1e3, np.real(R_tau_lag1))
+        for ind in range(1, proc_imgs.shape[0]):
+            plt.plot(az/1e3, np.real(proc_imgs[ind, n_lags+1, :]))
+
+    else:
+        proc_img = proc_imgs
+
+    #Estimate nrcs
+    nrcs_cosar = np.abs(proc_img[n_lags, :])
+    #FIXME
+    #Show results
+        #az in km
+    az = az / 1e3
+    #plt.figure()
+    #csl.plot_nrcs(az, nrcs, az, nrcs_cosar, dB=dB)
+    #Unwrap CoSAR using known true phase. Since this is not a phase unwrapping
+    #exercise, this is not cheating.
+    R_tau_lag0_CoSAR = proc_img[n_lags, :]
+    R_tau_lag0_CoSAR_a = np.abs(R_tau_lag0_CoSAR)
+    R_tau_lag0_CoSAR_p = (np.angle(R_tau_lag0_CoSAR * np.conj(R_tau_lag0)) +
+                          h_phase)
+    R_tau_lag1_CoSAR = proc_img[n_lags+1, :]
+    R_tau_lag1_CoSAR_a = np.abs(R_tau_lag1_CoSAR)
+    R_tau_lag1_CoSAR_p = (np.angle(R_tau_lag1_CoSAR * np.conj(R_tau_lag1)) +
+                          dop_phase + h_phase)
+    R_tau_lag1m_CoSAR = proc_img[n_lags-1, :]
+    R_tau_lag1m_CoSAR_a = np.abs(R_tau_lag1m_CoSAR)
+    R_tau_lag1m_CoSAR_p = (np.angle(R_tau_lag1m_CoSAR * np.conj(R_tau_lag1m)) -
+                           dop_phase + h_phase)
+
+    plt.figure()
+    csl.plot_R_tau(az, abs(R_tau_lag0), h_phase,
+                   az, R_tau_lag0_CoSAR_a, R_tau_lag0_CoSAR_p, dB=dB,
+                   title=r'$R_{\tau}(x,\tau=0)$', fontsize=fontsize)
+    if type(fileroot) == str:
+            filename = fileroot+'_R_tau0.eps'
+            plt.savefig(filename)
+            filename = fileroot+'_R_tau0.png'
+            plt.savefig(filename)
+
+    plt.figure()
+    csl.plot_R_tau(az, abs(R_tau_lag1), dop_phase + h_phase,
+                   az, R_tau_lag1_CoSAR_a, R_tau_lag1_CoSAR_p, dB=dB,
+                   title=r'$R_{\tau}(x,\tau=1/PRF)$', fontsize=fontsize)
+    if type(fileroot) == str:
+            filename = fileroot+'_R_tau1.eps'
+            plt.savefig(filename)
+            filename = fileroot+'_R_tau1.png'
+            plt.savefig(filename)
+
+    plt.figure()
+    csl.plot_R_tau(az, abs(R_tau_lag1m), -1.0 * dop_phase + h_phase,
+                   az, R_tau_lag1m_CoSAR_a, R_tau_lag1m_CoSAR_p, dB=dB,
+                   title=r'$R_{\tau}(x,\tau=-1/PRF)$', fontsize=fontsize)
+    if type(fileroot) == str:
+            filename = fileroot+'_R_tau-1.eps'
+            plt.savefig(filename)
+            filename = fileroot+'_R_tau-1.png'
+            plt.savefig(filename)
+
+
+def test_geosim1d():
+    """
+    Test routine for geosim1d
+    """
+    #Define geometry and radar
+    cs = csl.CoSARGeometry(0.8, 25, 3, 10e9, 10)
+    #Define ocean characteristics
+    sf = csl.CoSARSurface(0.1, 10, 1, 0, 1000)
+    #Run simulation
+    sim_res = geosim1d(cs, sf, 250.0, n_rg=100, PRF=10.0, Tint=300.0)
+
+

cosar/geosim1d_test.py

+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jan 10 10:29:26 2014
+
+@author: Paco López Dekker
+"""
+#%%
+from cosar.geosim1d import geosim1d_analyze, geosim1d, geosim1d_montecarlo,\
+                     geosim1d_save, geosim1d_load
+import cosar.simplesim as csl
+import os
+
+#%%
+n_it = 20
+v_stdev = 1.0  # Small to avoid big ambiguity issues
+h_stdev = 1.0
+dv = 6.0
+PRF = 200.0
+tau = 0.01
+Tint = 30.0
+
+respath = "/Users/plopezdekker/Documents/WORK/CoSAR/RESULTS/Figures"
+os.chdir(respath)
+#Define geometry and radar
+cs = csl.CoSARGeometry(0.8, 25, dv, 10e9, 10)
+#Define ocean characteristics
+sf = csl.CoSARSurface(tau, 10, v_stdev, h_stdev, 1000)
+#Run simulations
+Tint = 100
+sim_res = geosim1d_montecarlo(cs, sf, 250.0, n_rg=10, PRF=PRF, Tint=Tint,
+                              surfaceseed=1, n_it=n_it)
+geosim1d_save(sim_res, 'sim_PRF200_Tint100_nr10')
+sim_res = geosim1d_load('sim_PRF200_Tint100_nr10')
+geosim1d_analyze(cs, sf, sim_res, fileroot='sim_PRF200_Tint100_nr10', dB=True)
+
+Tint = 300
+#sim_res = geosim1d_montecarlo(cs, sf, 250.0, n_rg=10, PRF=PRF, Tint=Tint,
+#                              surfaceseed=1, n_it=n_it)
+#geosim1d_save(sim_res, 'sim_PRF200_Tint300_nr10')
+sim_res = geosim1d_load('sim_PRF200_Tint300_nr10')
+geosim1d_analyze(cs, sf, sim_res, fileroot='sim_PRF200_Tint300_nr10', dB=True)
+
+Tint = 600
+#sim_res = geosim1d_montecarlo(cs, sf, 250.0, n_rg=10, PRF=PRF, Tint=Tint,
+#                              surfaceseed=1, n_it=n_it)
+#geosim1d_save(sim_res, 'sim_PRF200_Tint600_nr10')
+sim_res = geosim1d_load('sim_PRF200_Tint600_nr10')
+geosim1d_analyze(cs, sf, sim_res, fileroot='sim_PRF200_Tint600_nr10', dB=True)
+
+# %%

cosar/results.py

+    # -*- coding: utf-8 -*-
+
+import sys
+import os
+import time 
+import numpy as np
+import cosar as cosar
+#from scipy import weave
+from scipy import signal
+import matplotlib.pyplot as plt
+
+
+def igarss2013(figs,respath=0):
+    """Generate figures for IGARSS"""
+    if type(respath) != str:
+        respath = "/Users/plopezdekker/Documents/WORK/CoSAR/RESULTS/Figures"
+    os.chdir(respath)
+    if figs == 1:
+        filename = "sar_coh_pt_64rg_T40.png"
+        cosar.simplesim(256,64,tau=-1,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename)
+        filename = "sar_tau20_pt_64rg_T40.png"
+        cosar.simplesim(256,64,tau=20,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename)
+        filename = "sar_tau10_pt_64rg_T40.png"
+        cosar.simplesim(256,64,tau=10,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename)
+        filename = "sar_tau5_pt_64rg_T40.png"
+        cosar.simplesim(256,64,tau=5,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename)
+        filename = "sar_tau1_pt_64rg_T40.png"
+        cosar.simplesim(256,64,tau=1,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename)
+        filename = "sar_tau01_pt_64rg_T40.png"
+        cosar.simplesim(256,64,tau=0.1,Tint=40.,dt=0,show_SAR=True,show_COSAR=False,dop_rms=0,pt=True,filesave=filename)
+    if figs == 2:
+        na = 256
+        nr = 256
+        Tint = 40.0
+        dt = 0
+        taus = [20,10,5,1,0.1]
+        filenames = ["sar_vs_cosar_tau20_pt_256rg_T40.png",
+                     "sar_vs_cosar_tau10_pt_256rg_T40.png",
+                     "sar_vs_cosar_tau5_pt_256rg_T40.png",
+                     "sar_vs_cosar_tau1_pt_256rg_T40.png",
+                     "sar_vs_cosar_tau01_pt_256rg_T40.png"]
+        for ind in range(0,len(taus)):
+            cosar.simplesim(na,nr,tau=taus[ind],Tint=Tint,dt=dt,show_SAR=True,show_COSAR=True,dop_rms=0,pt=True,filesave=filenames[ind])
+    if figs == 3:
+        filename = "sar_vs_cosar_tau01_seed1_1024rg_T20.png"
+        cosar.simplesim(256,1024,tau=0.01,Tint=20.,dt=0,show_SAR=True,show_COSAR=True,dop_rms=0.5,seed=1,filesave=filename)
+        nrs = [128,256,512,1024]
+        filenames = ["cosar_tau01_seed1_128rg_T20.png",
+                     "cosar_tau01_seed1_256rg_T20.png",
+                     "cosar_tau01_seed1_512rg_T20.png",
+                     "cosar_tau01_seed1_1024rg_T20.png"]
+        for ind in range(0,len(nrs)):
+            cosar.simplesim(256,nrs[ind],tau=0.01,Tint=20.,dt=0,show_SAR=False,show_COSAR=True,dop_rms=0.5,seed=1,filesave=filenames[ind])
+        
\ No newline at end of file

cosar/simplesim.py

+'''
+Created on May 31, 2013
+
+@author: lope_fr
+'''
+import sys
+import time
+import numpy as np
+#from scipy import weave
+from scipy import signal
+import matplotlib.pyplot as plt
+from scipy import constants
+from drama import utils
+from drama import geo as sar_geo
+#import IPython
+import matplotlib
+matplotlib.rcParams['ps.fonttype'] = 42
+#matplotlib.rcParams['text.usetex'] = True
+
+class CoSARGeometry:
+    """A class containing the relevant CoSAR orbital and geometricparameters"""
+    _orbit_height = 35786e3
+
+    def __init__(self, dkr, theta_i, dv, f0, antenna_length):
+        self.dkr = 1.0*dkr
+        self.theta_i = theta_i*np.pi/180
+        _dv = 1.0*np.array(dv)
+        if _dv.size == 1:
+            _dv = np.array([dv, 0.0, 0.0])
+        self.dv = _dv
+        self.f0 = 1.0*f0
+        self.L_antenna = 1.0 * antenna_length
+
+    def wavelength(self):
+        return constants.c/self.f0
+
+    def dtheta_i(self):
+        res = self.dkr * self.wavelength() / (4 * np.pi * np.cos(self.theta_i))
+        return res
+
+    def h_amb(self):
+        #res = (2 * np.pi * np.sin(self.theta_i) *
+        #       np.cos(self.theta_i) / self.dkr)
+        res = (2 * np.pi * np.cos(self.theta_i) /
+               (np.sin(self.theta_i) * self.dkr))
+        return res
+
+    def kz(self):
+        return 2 * np.pi / self.h_amb()
+
+    def R0(self):
+        return sar_geo.inc_to_sr(self.theta_i, self._orbit_height)
+
+    def L_illuminated(self):
+        return self.wavelength()/self.L_antenna * self.R0()
+
+    def dDoppler_bandwidth(self):
+        return 2*self.dv[0]/self.L_antenna
+
+
+class CoSARSurface:
+    """A class defining the scene"""
+    def __init__(self, tau, NRCS_stdev, vr_stdev, h_stdev,
+                 L_NRCS, L_vr=0, L_h=0,
+                 NRCS_mean=-10.0, vr_mean=0.0):
+        self.tau = tau
+        self.NRCS_stdev = NRCS_stdev
+        self.NRCS_mean = NRCS_mean
+        self.vr_stdev = vr_stdev
+        self.vr_mean = vr_stdev
+        self.h_stdev = h_stdev
+        self.L_NRCS = L_NRCS
+        if L_vr == 0:
+            self.L_vr = self.L_NRCS
+        else:
+            self.L_vr = L_vr
+        if L_h == 0:
+            self.L_h = 5.0 * self.L_NRCS
+        else:
+            self.L_h = L_h
+
+    def scene(self, L, dx, f0):
+        """Generate a realization of the scene"""
+        naz = int(np.round(L/dx))
+        az = np.linspace(-naz/2*dx, naz/2*dx, naz)
+        wavelength = constants.c/f0
+        dop_stdev = 2*self.vr_stdev/wavelength
+        dop_mean = 2*self.vr_mean/wavelength
+        nrcs = get_nrcs(naz, int(np.round(self.L_NRCS/dx)), pt=False,
+                        NRCS_mean=self.NRCS_mean, NRCS_stdev=self.NRCS_stdev)
+        dop = get_Doppler(naz, int(np.round(self.L_vr/dx)), dop_stdev)
+        dop = dop + dop_mean
+        #Reuse get_Doppler function to generate height profile
+        h = get_Doppler(naz, int(np.round(self.L_h/dx)), self.h_stdev)
+        nrcs.reshape(naz, 1)
+        dop.reshape(naz, 1)
+        h.reshape(naz, 1)
+        return (nrcs, dop, h, az)
+
+    def L_min(self):
+        """Returns the minimum length scale"""
+        return np.min(np.array([self.L_h, self.L_NRCS, self.L_vr]))
+
+
+def Q_w(R, NRCS, dx, N_r, T_i, tau_ca, az_res, F_n=1.0):
+    """
+    Estimate CoSAR imaging quality.
+    R is the space-dependent expected correlation
+    NRCS is the corresponding NRCS
+    dx is the grid spacing
+    N_r the number of range gates averaged
+    T_i the integration time
+    tau_ca the decorrelation time
+    az_res the azimuth resolution
+    """
+    Q_scaling = N_r / F_n * T_i / tau_ca * az_res**2
+    Power = np.sum(NRCS) * dx
+    Q_w = Q_scaling * np.abs(R)**2 / Power**2
+    return Q_w
+
+
+def get_nrcs(naz, gaussian_length, pt=False, SCR=30,
+             NRCS_mean=-10, NRCS_stdev=10.0):
+    """Generate a random but smooth NRCS"""
+    #We use a uniform distribution
+    dist_amp = np.sqrt(12) * NRCS_stdev
+    nrcs_db = dist_amp * np.random.random(naz) - dist_amp / 2 + NRCS_mean
+    #A Gaussian Filter
+    gfilt = signal.gaussian(int(naz / 4), gaussian_length, sym=True)
+    gfilt = gfilt / gfilt.sum()
+    nrcs_db_filt = signal.fftconvolve(nrcs_db, gfilt, mode='same')
+    nrcs_filt = utils.db2lin(nrcs_db_filt)
+    if pt:
+        nrcs_filt[:] = utils.db2lin(-SCR)
+        nrcs_filt[[naz / 4, naz / 2]] = 1
+    return nrcs_filt
+
+
+def get_Doppler(naz, gaussian_length, dop_rms=1):
+    """Generate random Doppler frequencies"""
+    dop = np.random.randn(naz)
+    #A Gaussian Filter
+    gfilt = signal.gaussian(int(naz / 4), gaussian_length, sym=True)
+    gfilt = gfilt / gfilt.sum()
+    dop_filt = signal.fftconvolve(dop, gfilt, mode='same')
+    norm_factor = dop_rms/np.std(dop_filt)
+    dop_filt = dop_filt * norm_factor
+    return dop_filt
+
+
+def cosar_proc_1d(rx1, rx2, pos1, pos2, az, k0, han=False, dt=0):
+    """CoSAR processing assuming no range migrations"""
+    #Cross correlate signals
+    if dt == 0:
+        corr_est = np.mean(rx2 * np.conjugate(rx1), axis=1)
+    else:
+        if dt > 0:
+            corr_p2p = rx2[dt:, :] * np.conjugate(rx1)[:-dt, :]
+        else:
+            corr_p2p = rx2[0:dt, :] * np.conjugate(rx1)[-dt:, :]
+        corr_est = np.mean(corr_p2p, axis=1)
+
+    naz = az.size
+    nt = corr_est.size
+    if han:
+        az_win = np.hanning(nt)
+    else:
+        az_win = np.ones(nt)
+
+    T_i_scaling = np.mean(az_win**2)
+    #Result
+    s_est = np.zeros(naz, dtype=np.complex)
+    for t_ind in range(0, corr_est.size):
+        #Calculate ranges from targets to both radars
+        #r_2 = ((pos1[t_ind, 1])**2 + (az - pos1[t_ind, 0])**2 +
+        #       pos1[t_ind, 2] ** 2 )
+        r1 = np.sqrt((pos1[t_ind, 1])**2 + (az - pos1[t_ind, 0])**2 +
+                     pos1[t_ind, 2]**2)
+
+        r2 = np.sqrt((pos2[t_ind + dt, 1])**2 +
+                     (az - pos2[t_ind + dt, 0])**2 +
+                     pos2[t_ind+dt, 2]**2)
+
+        dphase = 2*k0*(r2-r1)
+        phasor = np.exp(1j*dphase)
+        s_est = s_est + az_win[t_ind] * phasor * corr_est[t_ind]
+    #Normalize
+    s_est = s_est  # / corr_est.size / 2
+    return (s_est, T_i_scaling)
+
+
+def sar_proc_1d(rx, pos, az, k0, han=False):
+    """SAR processing assuming no range migrations using backprojection"""
+    nt = rx.shape[0]
+    nrg = rx.shape[1]
+    naz = az.size
+    if han:
+        az_win = np.hanning(nt)
+    else:
+        az_win = np.ones(nt)
+    #Results
+    T_i_scaling = np.mean(az_win**2)
+    s_est = np.zeros([naz, nrg], dtype=np.complex)
+    for t_ind in range(0, nt):
+        r = (np.sqrt((pos[t_ind, 1])**2 + (az - pos[t_ind, 0])**2 +
+             pos[t_ind, 2]**2))
+        phase = 2 * k0 * r
+        phasor = np.exp(1j * phase)
+        s_est = s_est + az_win[t_ind] * rx[t_ind, :] * phasor.reshape([naz, 1])
+    pow_est = np.mean(np.abs(s_est)**2, axis=1)
+    return (pow_est/nt, T_i_scaling)
+
+
+def plot_nrcs(az_true, nrcs_true, az_proc, nrcs_proc, dB=True, nrcs_err=0):
+    """Routine to plot cosar estimated NRCS and ground truth"""
+    ax = plt.subplot(1, 1, 1)
+    res_dim = nrcs_proc.shape
+    if np.size(res_dim) == 2:
+        n_res = res_dim[0]
+    else:
+        n_res = 1
+        nrcs_proc.shape = [n_res, res_dim[0]]
+    if dB:
+        p_true_nrcs, = plt.plot(az_true, utils.db(nrcs_true, norm=True), 'b',
+                                label="True")
+        p_cosar_nrcs, = plt.plot(az_proc, utils.db(nrcs_proc[0, :], norm=True),
+                                 'r', label="CoSAR")
+        for it_res in range(1, n_res):
+            plt.plot(az_proc, utils.db(nrcs_proc[it_res, :], norm=True), 'r')
+        if np.size(nrcs_err) == np.size(nrcs_true):
+            p_error, = plt.plot(az_true,
+                                utils.db(nrcs_true + nrcs_err, norm=True),
+                                'b--')
+
+        plt.xlabel('Azimuth [km]')
+        plt.ylabel('NRCS [dB]')
+    else:
+        p_true_nrcs, = plt.plot(az_true, nrcs_true, 'b',
+                                label="True")
+        p_cosar_nrcs, = plt.plot(az_proc, nrcs_proc[0, :], 'r',
+                                 label="CoSAR")
+        for it_res in range(1, n_res):
+            plt.plot(az_proc, nrcs_proc[it_res, :], 'r')
+        if np.size(nrcs_err) == np.size(nrcs_true):
+            p_error, = plt.plot(az_true,
+                                nrcs_true + nrcs_err, 'b--')
+        plt.xlabel('Azimuth [km]')
+        plt.ylabel('NRCS')
+    #Add labels
+    handles, labels = ax.get_legend_handles_labels()
+    plt.legend(handles, labels, loc=3)
+    plt.show()
+
+
+def plot_R_tau(az_true, r_true_a, r_true_p,
+               az_proc, r_proc_a, r_proc_p, dB=True, r_err=0,
+               title=False, fontsize=12):
+    """
+    Routine to plot CoSAR estimated R_tau and ground truth
+    Variables ending in _a are amplitude, those ending in _p are phases
+    """
+    ax = plt.subplot(2, 1, 1)
+    res_dim = r_proc_a.shape
+    if np.size(res_dim) == 2:
+        n_res = res_dim[0]
+    else:
+        n_res = 1
+        r_proc_a.shape = [n_res, res_dim[0]]
+        r_proc_p.shape = [n_res, res_dim[0]]
+    if dB:
+        p_true_r, = plt.plot(az_true, utils.db(r_true_a, norm=True), 'b',
+                             label="True")
+        p_cosar_r, = plt.plot(az_proc,
+                              utils.db(r_proc_a[0, :], norm=True),
+                              'r', label="CoSAR")
+        for it_res in range(1, n_res):
+            plt.plot(az_proc,
+                     utils.db(r_proc_a[it_res, :], norm=True),
+                     'r')
+        if np.size(r_err) == np.size(r_true_a):
+            p_error, = plt.plot(az_true,
+                                utils.db(r_true_a + r_err, norm=True),
+                                'b--')
+
+        #plt.xlabel('Azimuth [km]')
+        plt.ylabel(r'$|R|$ $\mathrm{[dB]}$', fontsize=fontsize)
+    else:
+        p_true_r, = plt.plot(az_true, r_true_a, 'b',
+                             label="True")
+        p_cosar_r, = plt.plot(az_proc, r_proc_a[0, :], 'r',
+                              label="CoSAR")
+        for it_res in range(1, n_res):
+            plt.plot(az_proc, r_proc_a[it_res, :], 'r')
+        if np.size(r_err) == np.size(r_true_a):
+            p_error, = plt.plot(az_true,
+                                r_true_a + r_err, 'b--')
+        #plt.xlabel('Azimuth [km]')
+        plt.ylabel(r'$|R|$', fontsize=fontsize)
+    if type(title) == str:
+        plt.title(title, fontsize=fontsize)
+    #Add labels
+    handles, labels = ax.get_legend_handles_labels()
+    plt.legend(handles, labels, loc=3, fontsize=fontsize)
+    #Now plot the phase
+    ax = plt.subplot(2, 1, 2)
+    p_true_r_ph, = plt.plot(az_true, r_true_p, 'b', label='True')
+    p_true_r_ph, = plt.plot(az_proc, r_proc_p[0, :], 'r', label='CoSAR')
+    for it_res in range(1, n_res):
+        plt.plot(az_proc, r_proc_p[it_res, :], 'r')
+
+    plt.xlabel('Azimuth [km]', fontsize=fontsize)
+    plt.ylabel(r'$\mathrm{arg}(R)$ $\mathrm{[rad]}$', fontsize=fontsize)
+
+
+
+def simplesim_show(az,true_nrcs,est_nrcs,true_Doppler=0,est_Doppler=0,sar_nrcs=0,
+                   show_SAR=False,show_Doppler=False,show_COSAR=True,filesave=0):
+    """Routine to plot simplesim simulation results"""
+    if show_Doppler:
+        ax = plt.subplot(2,1,1)
+    else:
+        ax = plt.subplot(1,1,1)
+    p_true_nrcs, = plt.plot(az,utils.db(true_nrcs,norm=True),'b',label="True")
+    if show_COSAR:
+        p_cosar_nrcs, = plt.plot(az,utils.db(est_nrcs,norm=True),'r',label="CoSAR")
+    plt.xlabel('Azimuth')
+    plt.ylabel('Intensity [dB]')
+    if show_SAR:
+        plt.plot(az,utils.db(sar_nrcs,norm=True),'k--',label="SAR")
+    #Add labels
+    handles, labels = ax.get_legend_handles_labels()
+    plt.legend(handles,labels,loc=3)
+    if show_Doppler > 0:
+        plt.subplot(2,1,2)
+        plt.plot(az,true_Doppler,'b',az,est_Doppler,'r')
+        plt.xlabel('Azimuth')
+        plt.ylabel('Doppler [Hz]')
+    plt.show()
+    #Check if we want to print this
+    if type(filesave) == str:
+        plt.savefig(filesave)
+        plt.close()
+
+def simplesim_scene_init(naz,length_scale,dop_rms=1.0,pt='False'):
+    nrcs = get_nrcs(naz,length_scale,pt=pt)
+    dop = get_Doppler(naz,10,dop_rms)
+    nrcs.reshape(naz,1)
+    dop.reshape(naz,1)
+    return (nrcs,dop)
+
+def simplesim(naz,nrg,tau=0.0,pt=False,hamming=True,dop_rms=1.0,Tint=20.0,dt=0,
+              show_SAR=False,show_COSAR=True,filesave=0,seed=None):
+
+    #Create some complexarray of random numbers
+
+    dx = 1.0
+    #Spacecraft height and ground range
+    z1 = z2 = 5e3
+    x1 = x2 = -5e3
+    R0 = 10e3
+    v1, v2 = 1.0, -1.0
+    f0 = 10e9
+    k0 = 2*np.pi*f0/constants.c
+    l0 = constants.c/f0
+    fDop_max = naz * dx / R0 * 2 * abs(v2-v1)/l0
+    PRF = 8*fDop_max
+    print("PRF = %i", int(PRF))
+    t_v = np.linspace(-Tint/2, Tint/2, Tint*PRF)
+    az = np.linspace(-naz/2*dx, naz/2*dx, naz)
+    ##Seed simulator if needed
+    np.random.seed(seed)
+    #nrcs
+    nrcs, dop = simplesim_scene_init(naz, 10, pt=pt)
+    dop_phasor = np.exp(2j*np.pi*dop/PRF)
+    #print dop_phasor.shape
+    expected_amp = np.sqrt(nrcs)
+    #teporal decorrelation factor or whatever
+    # Get trajectory of spacecraft
+    pos1 = np.zeros([Tint*PRF, 3])
+    pos2 = np.zeros([Tint*PRF, 3])
+    pos1[:, 0] = v1*t_v
+    pos2[:, 0] = v2*t_v
+    pos1[:, 1] = x1
+    pos2[:, 1] = x2
+    pos1[:, 2] = z1
+    pos2[:, 2] = z2
+    #Correlation coefficcient
+    # alpha**(tau*PRF) = 1/e
+    # tau*PRF * log(alpha)=-1 alpha = exp(-1/(tau*PRF)
+    alpha = 1
+    if tau > 0:
+        alpha = np.exp(-1.0/(tau*PRF))
+    if nrg > 1:
+        uscat = np.random.randn(naz, nrg) + 1j*np.random.randn(naz, nrg)
+    else:
+        uscat = np.exp(1j*2*np.pi*np.random.rand(naz, nrg))
+    for t_ind in range(0, t_v.size):
+        #Generate random complex scatterers
+        uscat = (alpha * uscat * dop_phasor.reshape(naz, 1) +
+                 np.sqrt(1-alpha**2) * (np.random.randn(naz, nrg) +
+                                        1j * np.random.randn(naz, nrg)))
+        #Multipliy them with a fixed amplitude
+        scat = uscat * expected_amp.reshape(naz, 1)
+        #Calculate ranges from targets to both radars
+        r1 = ((np.sqrt((pos1[t_ind, 1])**2 +
+              (az - pos1[t_ind, 0])**2 + pos1[t_ind, 2]**2)).reshape(naz, 1))
+        r2 = ((np.sqrt((pos2[t_ind, 1])**2 +
+              (az - pos2[t_ind, 0])**2 + pos2[t_ind, 2]**2)).reshape(naz, 1))
+        #Include 2-way phase in field
+        s1 = scat * np.exp(-2j*k0*r1)
+        s2 = scat * np.exp(-2j*k0*r2)
+        # Compute range profiles by integrating in azimuth.
+        # This assumed that there is no range migration
+        rx1_ = s1.sum(axis=0)
+        rx2_ = s2.sum(axis=0)
+        # Save in a temporal stack
+        if t_ind == 0:
+            rx1 = rx1_.copy()
+            rx2 = rx2_.copy()
+        else:
+            rx1 = np.vstack((rx1, rx1_))
+            rx2 = np.vstack((rx2, rx2_))
+    #Now we have a range-slow time stack of CoSAR data
+    proc_img, T_i_scaling = cosar_proc_1d(rx1, rx2, pos1, pos2, az, k0,
+                                          han=hamming)
+    if dt > 0:
+        proc_ximg = cosar_proc_1d(rx1, rx2, pos1, pos2, az, k0, han=hamming,
+                                  dt=dt)
+        #The phase of this should be the Doppler phase
+        est_Doppler = np.angle(proc_ximg)/(dt/PRF)/(2*np.pi)
+    else:
+        est_Doppler = 0
+
+    if show_SAR:
+        sar_img = sar_proc_1d(rx1, pos1, az, k0, han=hamming)
+    else:
+        sar_img = 0
+    #IPython.embed()
+    #Show results
+    #plt.ion()
+    simplesim_show(az, nrcs, np.abs(proc_img), dop, est_Doppler, sar_img,
+                   show_SAR=show_SAR, show_Doppler=(dt > 0),
+                   show_COSAR=show_COSAR, filesave=filesave)
+    print("Done!")
+

cosar/viciprep.py

+__author__ = "Paco Lopez Dekker"
+__email__ = "F.LopezDekker@tudeft.nl"
+
+import numpy as np
+import scipy.constants as cnst
+
+def spectral_shift(f0, theta_i, dtheta_i):
+    return f0 * dtheta_i / np.tan(theta_i)
+
+
+def CoSAR_airborne(target_res, h=2.5e3, Bhor=50, v_sar = 100, dv_sar=2, theta_i_deg=45, az_bmwdth_deg=30,
+                   f0=np.array([1.24e9, 5.4e9, 9.65e9]),
+                   Bw=np.array([200, 200, 500]) * 1e6,
+                   antenna_length=None):
+
+
+    wl0 = cnst.c / f0
+    theta_i = np.radians(theta_i_deg)
+    U = 5
+    tau_c = 3 * wl0 / U
+    Bn = np.cos(theta_i) * Bhor
+    R0 = h / np.cos(theta_i)
+    az_bmwdth = np.radians(az_bmwdth_deg)
+    if not antenna_length is None:
+        az_bmwdth = wl0 / antenna_length
+
+    # Spectral shift
+    dtheta_i = Bn / R0
+    df = spectral_shift(f0, theta_i, dtheta_i)
+    print("Spectral shift [MHz]")
+    print(df/1e6)
+    #az_bmwdth = np.radians(25)
+
+
+    width_gr = R0 *  np.sin(az_bmwdth)
+
+    Lap = wl0 * R0 / 2 / target_res
+    T_CoSAR_i = Lap / dv_sar
+
+    ## CoSAR Quality factor
+    Bw_common =  Bw - df
+    dgr = cnst.c / 2 / Bw_common / np.sin(theta_i)
+    Nr = target_res / dgr
+
+    Q = Nr * T_CoSAR_i / tau_c * target_res ** 2 / az_bmwdth
+
+    # Printing of the resuls
+    print("Slant range: %f" % R0)
+
+    np.set_printoptions(precision=3)
+    print("Azimuth antenna foot print:")
+    print(width_gr)
+    print("SAR focused res (decorrelation limited)")
+    print(wl0 / (v_sar * tau_c)* R0/2)
+    print("CoSAR aperture lengths")
+    print(Lap)
+    print("Ground Distance travelled")
+    print(T_CoSAR_i * v_sar)
+    print("CoSAR integation time")
+    print(T_CoSAR_i)
+    print("Number of range gates")
+    print(Nr.astype(np.int))
+    print("Quality factor (Q):")
+    print(Q)
+    print("Pseudo SCR (dB):")
+    print(5*np.log10(Q))
+    np.set_printoptions()
\ No newline at end of file

cosarsim.py

+#!/usr/bin/env python
+
+'''
+COSAR Simulator launcher
+
+Usage ::
+	./cosar (GUI version)
+	./cosar -nogui /path/to/parfile (Console version)
+	
+'''
+
+import os
+import sys
+import argparse
+import cosar
+
+# Check Python version
+if tuple(sys.version_info[:2]) < (2, 7):
+	raise Exception('Python %d.%d is not supported. Please use version 2.6 '
+                    'or greater.\n' % tuple(sys.version_info[:2]))
+
+def main(argv):
+	"""Launches COSAR Simulator"""
+	#Handle command line argument
+	parser = argparse.ArgumentParser(description="Simple CoSAR simulation")
+	parser.add_argument("-p", "--point_target", action="store_true")
+	parser.add_argument("-w", "--hamming", action="store_true")
+	parser.add_argument("-a", "--azimuth_points", type=int, default = 256)
+	parser.add_argument("-r", "--range_samples", type=int, default = 1024)
+	parser.add_argument("-c", "--correlation_coefficient", type=float,default = 0.0)
+	args = parser.parse_args()
+	naz = args.azimuth_points
+	nrg = args.range_samples
+	alpha = args.correlation_coefficient
+	cosar.simplesim(naz, nrg, alpha=alpha, pt=args.point_target, hamming=args.hamming)
+	
+	
+
+if __name__ == '__main__':
+	main(sys.argv)
\ No newline at end of file

setup.cfg

+[metadata]
+name = CoSAR
+version = 0.1
+author = Paco López-Dekker
+author_email = F.LopezDekker@tudelft.nl
+description = CoSAR simulation code
+long_description = file: README.rst
+long_description_content_type = text/x-rst
+# url = https://gitlab.tudelft.nl/plopezdekker/s1sea
+#project_urls =
+#    Bug Tracker = https://gitlab.tudelft.nl/plopezdekker/waddensar/-/issues
+classifiers =
+    Programming Language :: Python :: 3
+    License :: OSI Approved :: GNU General Public License v3 (GPLv3)
+    Operating System :: OS Independent
+    Intended Audience :: Science/Research
+    Topic :: Scientific/Engineering
+
+[options]
+packages = find:
+python_requires = >=3.9
+install_requires =
+        numpy>=1.20
+        scipy
+        numexpr>=2.7
+        matplotlib
+        drama

setup.py

+from setuptools import setup, find_packages
+
+setup(
+    name='cosar',
+    version='24.03.26',
+    # packages=['oceansar'],
+    packages=find_packages(exclude=['contrib', 'docs', 'tests*']),
+    # package_dir={'': 'oceansar'},
+    # url='url = https://gitlab.tudelft.nl/plopezdekker/s1sea',
+    license='GPL-3.0',
+    author='Paco Lopez Dekker',
+    author_email='F.LopezDekker@tudelft.nl',
+    description='',
+    install_requires=['numpy', 'scipy', 'matplotlib', 'drama']
+)