From 6f04edb47aa4b74782d46976e8b3635d2f9e3d51 Mon Sep 17 00:00:00 2001
From: knarf <knarf@1bdb13ef-5d69-4035-bb54-08abeb3aa7f1>
Date: Tue, 17 Jun 2014 16:10:59 +0000
Subject: add script to generate a mode spectum

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TOVSolver/trunk@153 1bdb13ef-5d69-4035-bb54-08abeb3aa7f1
---
 src/util/mode_spectrum.py | 73 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 73 insertions(+)
 create mode 100755 src/util/mode_spectrum.py

diff --git a/src/util/mode_spectrum.py b/src/util/mode_spectrum.py
new file mode 100755
index 0000000..4169841
--- /dev/null
+++ b/src/util/mode_spectrum.py
@@ -0,0 +1,73 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+# © Frank Löffler <knarf@cct.lsu.edu>
+# License: GPLv2+
+
+# Compute power spectral density of given data
+from plot_defaults import *
+from matplotlib.mlab import detrend_linear
+
+xlim = (0,10)
+ylim = (-140,-100)
+
+# load data
+Fx, Fy = np.loadtxt("hydrobase::rho.maximum.asc.bz2", comments="#", usecols=(1,2), unpack=True)
+dt = M_to_ms/(Fx[1]-Fx[0])
+
+# slice data if wanted
+start = 0.0
+end   = 1.
+Fx = Fx[int(start*len(Fx)): int(end*len(Fx))]
+Fy = Fy[int(start*len(Fy)): int(end*len(Fy))]
+
+# Plot basics
+fig = plt.figure(figsize=(6, 3))
+fig.subplots_adjust(top=0.92,bottom=0.16, left=0.11,right=0.98)
+ax = fig.add_subplot(1,1,1)
+
+# mode names and frequencies for non-rotating Gamma=2 K=100 TOV stars
+modes = {
+"F" : 1.44252691528028,
+"H1": 3.95408396916149,
+"H2": 5.91494894170517,
+"H3": 7.77382493405710,
+"H4": 9.58768293806276,
+"H5": 11.3772129983097,
+"H6": 13.1520241905666,
+"H7": 14.9172321735655,
+}
+
+# plot modes as vertical lines and label them
+for mode, freq in modes.items():
+  ax.plot((freq,freq), ylim, 'b--')
+  ax.text(freq-0.1, ylim[1]+1, mode, fontsize=fontsize)
+#ax.plot((2*modes["F"], 2*modes["F"]), (-200,0), 'g--')
+#ax.plot((3*modes["F"], 3*modes["F"]), (-200,0), 'g--')
+#ax.plot((4*modes["F"], 4*modes["F"]), (-200,0), 'g--')
+
+# plot PSD
+frac = 0.5 * len(Fy)
+ax.psd(Fy, NFFT=int(frac), pad_to=5*int(frac), noverlap=int(frac*0.5), detrend=detrend_linear, Fs=dt, linestyle='-', color='black', scale_by_freq=False)
+#ax.psd(Fy, NFFT=int(len(Fy)*frac), pad_to=10*int(len(Fy)*frac), noverlap=int(len(Fy)*frac*0.50), detrend=detrend_linear, Fs=M_to_ms/(Fx[1]-Fx[0]), linestyle='-', marker='x', color='red')
+
+# plot properties
+ax.set_xlim(xlim)
+ax.set_ylim(ylim) 
+
+#plt.title('PSD', fontsize=fontsize)
+ax.set_xlabel(r'$f$ [kHz]')
+ax.set_ylabel(r'PSD [dB/Hz]')
+
+ax.xaxis.set_major_locator(mticker.MaxNLocator(10))
+ax.xaxis.set_minor_locator(mticker.MaxNLocator(20))
+ax.xaxis.grid(False)
+ax.yaxis.set_major_locator(mticker.MaxNLocator(5))
+ax.yaxis.set_minor_locator(mticker.MaxNLocator(9))
+ax.yaxis.grid(False)
+set_tick_sizes(ax, 8, 4)
+
+
+#plt.show()
+plt.savefig('mode_spectrum.pdf')
+
-- 
cgit v1.2.3