diff --git a/AltitudeProfile.py b/AltitudeProfile.py
index 7e976d6..11bfd3c 100755
--- a/AltitudeProfile.py
+++ b/AltitudeProfile.py
@@ -1,38 +1,17 @@
 #!/usr/bin/env python
 """ Height Profile Example """
-import pyiri2016
+import pyiri2016 as iri
+import pyiri2016.plots as piri
 #
 import numpy as np
-from matplotlib.pyplot import figure, show
+from matplotlib.pyplot import  show
 
 glat, glon = -11.95, -76.77
 
 alt_km = np.arange(80,1000,20.)
 
-iri = pyiri2016.IRI('2012-08-21T12', alt_km, glat, glon)
+iono = iri.IRI('2012-08-21T12', alt_km, glat, glon)
 
-fig = figure(figsize=(16,6))
-axs = fig.subplots(1,2)
-
-fig.suptitle(f'{str(iri.time[0].values)[:-13]}\n Glat, Glon: {iri.glat}, {iri.attrs["glon"]}')
-
-
-pn = axs[0]
-pn.plot(iri['ne'].squeeze(), iri.alt_km, label='N$_e$')
-#pn.set_title(iri2016Obj.title1)
-pn.set_xlabel('Density (m$^{-3}$)')
-pn.set_ylabel('Altitude (km)')
-pn.set_xscale('log')
-pn.legend(loc='best')
-pn.grid(True)
-
-pn = axs[1]
-pn.plot(iri['Ti'].squeeze(), iri.alt_km, label='T$_i$')
-pn.plot(iri['Te'].squeeze(), iri.alt_km, label='T$_e$')
-#pn.set_title(iri2016Obj.title2)
-pn.set_xlabel('Temperature (K)')
-pn.set_ylabel('Altitude (km)')
-pn.legend(loc='best')
-pn.grid(True)
+piri.altprofile(iono)
 
 show()
\ No newline at end of file
diff --git a/LatitudeProfile.py b/LatitudeProfile.py
index 3ef3fe8..a4f48fa 100755
--- a/LatitudeProfile.py
+++ b/LatitudeProfile.py
@@ -1,45 +1,23 @@
 #!/usr/bin/env python
-import pyiri2016
-from numpy import arange
-from matplotlib.pyplot import figure,  show
+import pyiri2016 as iri
+import pyiri2016.plots as piri
+from matplotlib.pyplot import show
 
 """ Geog. Latitude Profile Example """
 
-latlim = [-60, 60]
-latstp = 2.
-sim = pyiri2016.geoprofile(altkm=600, latlim=latlim, dlat=latstp, \
-glon=-76.77,   time='2004-01-01T17')
 
-latbins = arange(latlim[0], latlim[1], latstp)
+if __name__ == '__main__':
+    from argparse import ArgumentParser
+    p = ArgumentParser()
+    p.add_argument('-o','--outfn',help='write data to file')
+    p = p.parse_args()
 
+    latlim = [-60, 60]
+    latstp = 2.
+    iono = iri.geoprofile(altkm=600, latlim=latlim, dlat=latstp,
+                          glon=-76.77,   time='2004-01-01T17')
 
-fig = figure(figsize=(8,12))
-axs = fig.subplots(2,1, sharex=True)
+    piri.latprofile(iono)
 
-pn = axs[0]
 
-pn.plot(latbins, sim['NmF2'].squeeze(), label='N$_m$F$_2$')
-pn.plot(latbins, sim['NmF1'].squeeze(), label='N$_m$F$_1$')
-pn.plot(latbins, sim['NmE'].squeeze(), label='N$_m$E')
-pn.set_title(str(sim.time[0].values)[:-13] + '  latitude'+str(latlim))
-pn.set_xlim(latbins[[0, -1]])
-pn.set_xlabel('Geog. Lat. ($^\circ$)')
-pn.set_ylabel('(m$^{-3}$)')
-pn.set_yscale('log')
-
-
-pn = axs[1]
-pn.plot(latbins, sim['hmF2'].squeeze(), label='h$_m$F$_2$')
-pn.plot(latbins, sim['hmF1'].squeeze(), label='h$_m$F$_1$')
-pn.plot(latbins, sim['hmE'].squeeze(), label='h$_m$E')
-pn.set_xlim(latbins[[0, -1]])
-pn.set_title(str(sim.time[0].values)[:-13] + '  latitude'+str(latlim))
-pn.set_xlabel('Geog. Lat. ($^\circ$)')
-pn.set_ylabel('(km)')
-
-for a in axs:
-    a.legend(loc='best')
-    a.grid(True)
-
-
-show()
+    show()
diff --git a/TimeProfile.py b/TimeProfile.py
index 26729a9..608871a 100755
--- a/TimeProfile.py
+++ b/TimeProfile.py
@@ -2,13 +2,10 @@
 """ Time Profile: IRI2016 """
 import numpy as np
 from datetime import timedelta
-from matplotlib.pyplot import figure, show
-try:
-    import ephem
-except ImportError:
-    ephem = None
+from matplotlib.pyplot import show
 #
-import pyiri2016
+import pyiri2016 as iri
+import pyiri2016.plots as piri
 
 
 # %% user parameters
@@ -17,88 +14,9 @@
 glat,glon = 65,-148
 alt_km = np.arange(120, 180, 20)
 # %% ru
-sim = pyiri2016.timeprofile(('2012-08-21','2012-08-22'),timedelta(hours=0.25),
+sim = iri.timeprofile(('2012-08-21','2012-08-22'),timedelta(hours=0.25),
                             alt_km,glat,glon)
 
-# %% Plots
-Nplot=3
-
-if Nplot>2:
-    fig = figure(figsize=(16,12))
-    axs = fig.subplots(3,1, sharex=True).ravel()
-else:
-    fig = figure(figsize=(16,6))
-    axs = fig.subplots(1,2).ravel()
-
-fig.suptitle(f'{str(sim.time[0].values)[:-13]} to {str(sim.time[-1].values)[:-13]}\n Glat, Glon: {sim.glat}, {sim.glon}')
-
-ax = axs[0]
-#NmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[2, :])
-#NmE = sim.loc[4, :]
-ax.plot(sim.time, sim['NmF2'].squeeze(), label='N$_m$F$_2$')
-#ax.plot(sim.time, NmF1, label='N$_m$F$_1$')
-#ax.plot(sim.time, NmE, label='N$_m$E')
-ax.set_title('Maximum number densities vs. ionospheric layer')
-ax.set_xlabel('Hour (UT)')
-ax.set_ylabel('(m$^{-3}$)')
-ax.set_yscale('log')
-ax.legend(loc='best')
-
-ax = axs[1]
-#hmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[3, :])
-#hmE = sim.b[5, :]
-ax.plot(sim.time, sim['hmF2'].squeeze(), label='h$_m$F$_2$')
-#ax.plot(sim.time, hmF1, label='h$_m$F$_1$')
-#ax.plot(sim.time, hmE, label='h$_m$E')
-ax.set_title('Height of maximum density vs. ionospheric layer')
-ax.set_xlabel('Hour (UT)')
-ax.set_ylabel('(km)')
-ax.legend(loc='best')
-# %%
-if Nplot > 2:
-    ax = axs[2]
-
-    for a in sim.alt_km:
-        ax.plot(sim.time, sim['ne'].squeeze(), marker='.', label=f'{a.item()} km')
-    ax.set_xlabel('time UTC (hours)')
-    ax.set_ylabel('[m$^{-3}$]')
-    ax.set_title(f'$N_e$ vs. altitude and time')
-    ax.set_yscale('log')
-    ax.legend(loc='best')
-# %%
-if Nplot > 4:
-    ax = axs[4]
-    tec = sim.b[36, :]
-    ax.plot(sim.time, tec, label=r'TEC')
-    ax.set_xlabel('Hour (UT)')
-    ax.set_ylabel('(m$^{-2}$)')
-    #ax.set_yscale('log')
-    ax.legend(loc='best')
-
-    ax = axs[5]
-    vy = sim.b[43, :]
-    ax.plot(sim.time, vy, label=r'V$_y$')
-    ax.set_xlabel('Hour (UT)')
-    ax.set_ylabel('(m/s)')
-    ax.legend(loc='best')
-
-for a in axs.ravel():
-    a.grid(True)
-
-
-if __name__ == '__main__':
-    from argparse import ArgumentParser
-    p = ArgumentParser(description='IRI2016 time profile plot')
-    p.add_argument('-t','--trange',help='START STOP STEP (hours) time [UTC]',nargs=3,
-                   default=('2012-08-21','2012-08-22',0.25))
-    p.add_argument('--alt',help='START STOP STEP altitude [km]',type=float, nargs=3,default=(120,180,20))
-    p.add_argument('-c','--latlon',help='geodetic coordinates of simulation',
-                   type=float,default=(65,-147.5))
-    p.add_argument('--f107',type=float,default=200.)
-    p.add_argument('--f107a', type=float,default=200.)
-    p.add_argument('--ap', type=int, default=4)
-    p.add_argument('--species',help='species to plot',nargs='+',default=('ne'))
-    p = p.parse_args()
-
+piri.timeprofile(sim)
 
 show()
diff --git a/pyiri2016/plots.py b/pyiri2016/plots.py
new file mode 100644
index 0000000..9160084
--- /dev/null
+++ b/pyiri2016/plots.py
@@ -0,0 +1,121 @@
+import xarray
+from matplotlib.pyplot import figure
+
+def timeprofile(iono:xarray.Dataset):
+    # %% Plots
+    Nplot=3
+
+    if Nplot>2:
+        fig = figure(figsize=(16,12))
+        axs = fig.subplots(3,1, sharex=True).ravel()
+    else:
+        fig = figure(figsize=(16,6))
+        axs = fig.subplots(1,2).ravel()
+
+    fig.suptitle(f'{str(iono.time[0].values)[:-13]} to {str(iono.time[-1].values)[:-13]}\n Glat, Glon: {iono.glat}, {iono.glon}')
+
+    ax = axs[0]
+
+    ax.plot(iono.time, iono['NmF2'].squeeze(), label='N$_m$F$_2$')
+    ax.plot(iono.time, iono['NmF1'].squeeze(), label='N$_m$F$_1$')
+    ax.plot(iono.time, iono['NmE'].squeeze(), label='N$_m$E')
+    ax.set_title('Maximum number densities vs. ionospheric layer')
+    ax.set_xlabel('Hour (UT)')
+    ax.set_ylabel('(m$^{-3}$)')
+    ax.set_yscale('log')
+    ax.legend(loc='best')
+
+    ax = axs[1]
+    ax.plot(iono.time, iono['hmF2'].squeeze(), label='h$_m$F$_2$')
+    ax.plot(iono.time, iono['hmF1'].squeeze(), label='h$_m$F$_1$')
+    ax.plot(iono.time, iono['hmE'].squeeze(), label='h$_m$E')
+    ax.set_title('Height of maximum density vs. ionospheric layer')
+    ax.set_xlabel('Hour (UT)')
+    ax.set_ylabel('(km)')
+    ax.legend(loc='best')
+    # %%
+    if Nplot > 2:
+        ax = axs[2]
+
+        for a in iono.alt_km:
+            ax.plot(iono.time, iono['ne'].squeeze(), marker='.', label=f'{a.item()} km')
+        ax.set_xlabel('time UTC (hours)')
+        ax.set_ylabel('[m$^{-3}$]')
+        ax.set_title(f'$N_e$ vs. altitude and time')
+        ax.set_yscale('log')
+        ax.legend(loc='best')
+    # %%
+    if Nplot > 4:
+        ax = axs[4]
+        tec = iono.b[36, :]
+        ax.plot(iono.time, tec, label=r'TEC')
+        ax.set_xlabel('Hour (UT)')
+        ax.set_ylabel('(m$^{-2}$)')
+        #ax.set_yscale('log')
+        ax.legend(loc='best')
+
+        ax = axs[5]
+        vy = iono.b[43, :]
+        ax.plot(iono.time, vy, label=r'V$_y$')
+        ax.set_xlabel('Hour (UT)')
+        ax.set_ylabel('(m/s)')
+        ax.legend(loc='best')
+
+    for a in axs.ravel():
+        a.grid(True)
+
+
+def altprofile(iono:xarray.Dataset):
+    fig = figure(figsize=(16,6))
+    axs = fig.subplots(1,2)
+
+    fig.suptitle(f'{str(iono.time[0].values)[:-13]}\n Glat, Glon: {iono.glat}, {iono.attrs["glon"]}')
+
+
+    pn = axs[0]
+    pn.plot(iono['ne'].squeeze(), iono.alt_km, label='N$_e$')
+    #pn.set_title(iri2016Obj.title1)
+    pn.set_xlabel('Density (m$^{-3}$)')
+    pn.set_ylabel('Altitude (km)')
+    pn.set_xscale('log')
+    pn.legend(loc='best')
+    pn.grid(True)
+
+    pn = axs[1]
+    pn.plot(iono['Ti'].squeeze(), iono.alt_km, label='T$_i$')
+    pn.plot(iono['Te'].squeeze(), iono.alt_km, label='T$_e$')
+    #pn.set_title(iri2016Obj.title2)
+    pn.set_xlabel('Temperature (K)')
+    pn.set_ylabel('Altitude (km)')
+    pn.legend(loc='best')
+    pn.grid(True)
+
+def latprofile(iono:xarray.Dataset):
+
+    fig = figure(figsize=(8,12))
+    axs = fig.subplots(2,1, sharex=True)
+
+    ax = axs[0]
+
+    ax.plot(iono.lat, iono['NmF2'].squeeze(), label='N$_m$F$_2$')
+    ax.plot(iono.lat, iono['NmF1'].squeeze(), label='N$_m$F$_1$')
+    ax.plot(iono.lat, iono['NmE'].squeeze(), label='N$_m$E')
+    ax.set_title(str(iono.time[0].values)[:-13] + f'  latitude {iono.lat[[0, -1]].values}')
+    #ax.set_xlim(iono.lat[[0, -1]])
+    ax.set_xlabel('Geog. Lat. ($^\circ$)')
+    ax.set_ylabel('(m$^{-3}$)')
+    ax.set_yscale('log')
+
+
+    ax = axs[1]
+    ax.plot(iono.lat, iono['hmF2'].squeeze(), label='h$_m$F$_2$')
+    ax.plot(iono.lat, iono['hmF1'].squeeze(), label='h$_m$F$_1$')
+    ax.plot(iono.lat, iono['hmE'].squeeze(), label='h$_m$E')
+    ax.set_xlim(iono.lat[[0, -1]])
+    ax.set_title(str(iono.time[0].values)[:-13] + f'  latitude  {iono.lat[[0, -1]].values}')
+    ax.set_xlabel('Geog. Lat. ($^\circ$)')
+    ax.set_ylabel('(km)')
+
+    for a in axs:
+        a.legend(loc='best')
+        a.grid(True)