User: knarf Date: 2011/09/12 11:29 PM
Added: /examples/tov/ mode_spectrum.py
Removed: /examples/tov/ mode_labels.sh, mode_spectrum.plot, mode_spectrum.svg, modes.dat
Modified: /examples/tov/ make_plots.sh, mode_spectrum.pdf
Log: use matplotlib for mode plot, calculate psd on the fly, including overlap with Hanning window, padding and linear detrending
Directory Changes:
Directory: /svn:executable/ ===========================
+ *
File Changes:
Directory: /examples/tov/ =========================
File [modified]: make_plots.sh Delta lines: +1 -2 =================================================================== --- examples/tov/make_plots.sh 2011-09-09 18:23:23 UTC (rev 130) +++ examples/tov/make_plots.sh 2011-09-13 04:29:02 UTC (rev 131) @@ -1,4 +1,3 @@ -gnuplot < mode_spectrum.plot -inkscape -A mode_spectrum.pdf mode_spectrum.svg +./mode_spectrum.py gnuplot < rho_max.plot inkscape -A rho_max.pdf rho_max.svg
File [removed]: mode_labels.sh Delta lines: +0 -1 =================================================================== --- examples/tov/mode_labels.sh 2011-09-09 18:23:23 UTC (rev 130) +++ examples/tov/mode_labels.sh 2011-09-13 04:29:02 UTC (rev 131) @@ -1 +0,0 @@ -perl -ane 'if ($F[2]) { print "set label "$F[2]" at ".($F[0]-0.15).",2.e-6\n"; }' modes.dat
File [modified]: mode_spectrum.pdf Delta lines: +0 -0 =================================================================== (Binary files differ)
File [removed]: mode_spectrum.plot Delta lines: +0 -13 =================================================================== --- examples/tov/mode_spectrum.plot 2011-09-09 18:23:23 UTC (rev 130) +++ examples/tov/mode_spectrum.plot 2011-09-13 04:29:02 UTC (rev 131) @@ -1,13 +0,0 @@ -set size ratio 0.3 -set term svg enhanced size 600,250 -set output 'mode_spectrum.svg' -set xlabel 'f/[kHz]' -set xtics 2 nomirror -set mxtics 2 -set x2tics scale 0 -set ylabel 'PSD' -set ytics 10 -#load '<./mode_labels.sh' -plot [:10][:-50] '<bzcat tov_2/hydrobase::rho.maximum.asc.spectrum.bz2' u ($1/1000):(10*log10($2)) w lp t '', \ -'<sed -e ''s/(.*)/\1 -80\n\1 -50\n/g'' modes.dat' u 1:3 w l t '', \ -'modes.dat' u 1:(0):x2ticlabels(2) w l t ''
File [added]: mode_spectrum.py Delta lines: +105 -0 =================================================================== --- examples/tov/mode_spectrum.py (rev 0) +++ examples/tov/mode_spectrum.py 2011-09-13 04:29:02 UTC (rev 131) @@ -0,0 +1,105 @@ +#!/usr/bin/python + +# Compute power spectral density of given data +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.ticker as mticker + +from matplotlib import rc +from matplotlib.mlab import detrend_linear + +# stuff +fontsize = 20 +xlim = (0,10) +ylim = (-140,-100) +rc('text', usetex=True) +rc('font', family='serif') +rc('font', serif='palatino') +rc('font', weight='bolder') +rc('mathtext', default='sf') +rc("lines", markeredgewidth=1) +rc("lines", linewidth=3) +rc('axes', labelsize=fontsize) +rc("axes", linewidth=2) +rc('xtick', labelsize=fontsize) +rc('ytick', labelsize=fontsize) +rc('legend', fontsize=fontsize) +rc('xtick.major', pad=8) +rc('ytick.major', pad=8) + +def set_tick_sizes(ax, major, minor): + for l in ax.get_xticklines() + ax.get_yticklines(): + l.set_markersize(major) + for tick in ax.xaxis.get_minor_ticks() + ax.yaxis.get_minor_ticks(): + tick.tick1line.set_markersize(minor) + tick.tick2line.set_markersize(minor) + ax.xaxis.LABELPAD=10. + ax.xaxis.OFFSETTEXTPAD=10. + +# constants +G = 6.673e-11 +c = 299792458 +M_sol = 1.98892e30 +# convertion factors +M_to_ms = 1./(1000*M_sol*G/(c*c*c)) + +# load data +Fx, Fy = np.loadtxt("tov_2/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() +fig.subplots_adjust(hspace=0.45, wspace=0.3, left=0.15, bottom=0.15) +ax = fig.add_subplot(1,1,1) + +# mode names and frequencies +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') +
Property changes on: examples/tov/mode_spectrum.py ___________________________________________________________________
File [removed]: mode_spectrum.svg Delta lines: +0 -210 =================================================================== --- examples/tov/mode_spectrum.svg 2011-09-09 18:23:23 UTC (rev 130) +++ examples/tov/mode_spectrum.svg 2011-09-13 04:29:02 UTC (rev 131) @@ -1,210 +0,0 @@ -<?xml version="1.0" encoding="utf-8" standalone="no"?> -<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" - "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd"> -<svg width="600" height="250" viewBox="0 0 600 250" - xmlns="http://www.w3.org/2000/svg" - xmlns:xlink="http://www.w3.org/1999/xlink%22%3E -<desc>Produced by GNUPLOT 4.2 patchlevel 4 </desc> - -<defs> - - - <circle id='gpDot' r='1'/> - <path id='gpPt0' style='stroke-width:0.222' d='M-1,0 h2 M0,-1 v2'/> - <path id='gpPt1' style='stroke-width:0.222' d='M-1,-1 L1,1 M1,-1 L-1,1'/> - <path id='gpPt2' style='stroke-width:0.222' d='M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1'/> - <rect id='gpPt3' style='stroke-width:0.222' x='-1' y='-1' width='2' height='2'/> - <use xlink:href='#gpPt3' id='gpPt4' style='fill:currentColor; stroke:none'/> - <circle id='gpPt5' style='stroke-width:0.222' cx='0' cy='0' r='1'/> - <use xlink:href='#gpPt5' id='gpPt6' style='fill:currentColor; stroke:none'/> - <path id='gpPt7' style='stroke-width:0.222' d='M0,-1.33 L-1.33,0.67 L1.33,0.67 z'/> - <use xlink:href='#gpPt7' id='gpPt8' style='fill:currentColor; stroke:none'/> - <use xlink:href='#gpPt7' id='gpPt9' transform='rotate(180)'/> - <use xlink:href='#gpPt9' id='gpPt10' style='fill:currentColor; stroke:none'/> - <use xlink:href='#gpPt3' id='gpPt11' transform='rotate(45)'/> - <use xlink:href='#gpPt11' id='gpPt12' style='fill:currentColor; stroke:none'/> -</defs> -<g style="fill:none; color:white; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter"> -</g> -<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter"> - <path d='M93.4,187.6 L102.4,187.6 M582.7,187.6 L573.7,187.6 '></path> - <g transform="translate(85.1,192.1)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:end"> - <text><tspan>-100</tspan> - </text> - </g> - <path d='M93.4,158.2 L102.4,158.2 M582.7,158.2 L573.7,158.2 '></path> - <g transform="translate(85.1,162.7)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:end"> - <text><tspan>-90</tspan> - </text> - </g> - <path d='M93.4,128.9 L102.4,128.9 M582.7,128.9 L573.7,128.9 '></path> - <g transform="translate(85.1,133.4)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:end"> - <text><tspan>-80</tspan> - </text> - </g> - <path d='M93.4,99.5 L102.4,99.5 M582.7,99.5 L573.7,99.5 '></path> - <g transform="translate(85.1,104.0)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:end"> - <text><tspan>-70</tspan> - </text> - </g> - <path d='M93.4,70.2 L102.4,70.2 M582.7,70.2 L573.7,70.2 '></path> - <g transform="translate(85.1,74.7)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:end"> - <text><tspan>-60</tspan> - </text> - </g> - <path d='M93.4,40.8 L102.4,40.8 M582.7,40.8 L573.7,40.8 '></path> - <g transform="translate(85.1,45.3)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:end"> - <text><tspan>-50</tspan> - </text> - </g> - <path d='M93.4,187.6 L93.4,178.6 '></path> - <g transform="translate(93.4,210.1)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan> 0</tspan> - </text> - </g> - <path d='M142.3,187.6 L142.3,183.1 M191.3,187.6 L191.3,178.6 '></path> - <g transform="translate(191.3,210.1)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan> 2</tspan> - </text> - </g> - <path d='M240.2,187.6 L240.2,183.1 M289.1,187.6 L289.1,178.6 '></path> - <g transform="translate(289.1,210.1)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan> 4</tspan> - </text> - </g> - <path d='M338.1,187.6 L338.1,183.1 M387.0,187.6 L387.0,178.6 '></path> - <g transform="translate(387.0,210.1)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan> 6</tspan> - </text> - </g> - <path d='M435.9,187.6 L435.9,183.1 M484.8,187.6 L484.8,178.6 '></path> - <g transform="translate(484.8,210.1)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan> 8</tspan> - </text> - </g> - <path d='M533.8,187.6 L533.8,183.1 M582.7,187.6 L582.7,178.6 '></path> - <g transform="translate(582.7,210.1)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan> 10</tspan> - </text> - </g> - <path d='M164.0,40.8 '></path> - <g transform="translate(164.0,27.3)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan>F</tspan> - </text> - </g> - <path d='M286.9,40.8 '></path> - <g transform="translate(286.9,27.3)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan>H1</tspan> - </text> - </g> - <path d='M382.8,40.8 '></path> - <g transform="translate(382.8,27.3)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan>H2</tspan> - </text> - </g> - <path d='M473.8,40.8 '></path> - <g transform="translate(473.8,27.3)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan>H3</tspan> - </text> - </g> - <path d='M562.5,40.8 '></path> - <g transform="translate(562.5,27.3)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan>H4</tspan> - </text> - </g> - <path d='M93.4,40.8 L93.4,187.6 L582.7,187.6 L582.7,40.8 L93.4,40.8 Z '></path> - <g transform="translate(39.1,114.2) rotate(-90)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan>PSD</tspan> - </text> - </g> - <g transform="translate(338.0,237.1)" style="stroke:none; fill:black; font-family:Arial; font-size:12.00pt; text-anchor:middle"> - <text><tspan>f/[kHz]</tspan> - </text> - </g> -</g> -<g style="fill:none; color:red; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter"> - <path d='M100.3,40.8 L101.0,63.7 L108.6,98.6 L116.2,106.9 L123.8,116.6 L131.4,118.4 L138.9,128.3 L146.5,122.9 - L154.1,101.9 L161.7,72.9 L169.3,102.9 L176.9,132.7 L184.5,137.7 L192.1,142.0 L199.7,152.4 L207.3,177.6 - L214.8,156.7 L222.4,187.4 L230.0,146.2 L237.6,133.1 L245.2,132.1 L252.8,135.5 L260.4,136.4 L268.0,119.8 - L275.6,106.7 L283.2,84.9 L290.7,83.4 L298.3,112.3 L305.9,113.1 L313.5,134.8 L321.1,121.9 L328.7,145.0 - L336.3,147.1 L343.9,134.9 L351.5,124.0 L359.1,117.5 L366.6,119.9 L374.2,100.8 L381.8,99.3 L389.4,96.2 - L397.0,130.7 L404.6,130.2 L412.2,130.3 L419.8,120.0 L427.4,148.9 L435.0,125.2 L442.6,129.9 L450.1,135.5 - L457.7,107.5 L465.3,116.5 L472.9,96.1 L480.5,119.5 L488.1,108.2 L495.7,111.9 L503.3,116.3 L510.9,122.5 - L518.5,122.2 L526.0,119.9 L533.6,127.1 L541.2,121.8 L548.8,125.6 L556.4,112.5 L564.0,118.2 L571.6,113.8 - L579.2,116.6 L582.7,117.5 '></path> - <use xlink:href='#gpPt0' transform='translate(101.0,63.7) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(108.6,98.6) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(116.2,106.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(123.8,116.6) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(131.4,118.4) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(138.9,128.3) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(146.5,122.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(154.1,101.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(161.7,72.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(169.3,102.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(176.9,132.7) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(184.5,137.7) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(192.1,142.0) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(199.7,152.4) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(207.3,177.6) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(214.8,156.7) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(222.4,187.4) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(230.0,146.2) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(237.6,133.1) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(245.2,132.1) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(252.8,135.5) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(260.4,136.4) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(268.0,119.8) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(275.6,106.7) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(283.2,84.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(290.7,83.4) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(298.3,112.3) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(305.9,113.1) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(313.5,134.8) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(321.1,121.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(328.7,145.0) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(336.3,147.1) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(343.9,134.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(351.5,124.0) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(359.1,117.5) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(366.6,119.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(374.2,100.8) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(381.8,99.3) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(389.4,96.2) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(397.0,130.7) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(404.6,130.2) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(412.2,130.3) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(419.8,120.0) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(427.4,148.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(435.0,125.2) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(442.6,129.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(450.1,135.5) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(457.7,107.5) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(465.3,116.5) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(472.9,96.1) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(480.5,119.5) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(488.1,108.2) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(495.7,111.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(503.3,116.3) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(510.9,122.5) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(518.5,122.2) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(526.0,119.9) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(533.6,127.1) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(541.2,121.8) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(548.8,125.6) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(556.4,112.5) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(564.0,118.2) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(571.6,113.8) scale(4.50)'/> - <use xlink:href='#gpPt0' transform='translate(579.2,116.6) scale(4.50)'/> -</g> -<g style="fill:none; color:green; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter"> - <path d='M164.0,128.9 L164.0,40.8 M286.9,128.9 L286.9,40.8 M382.8,128.9 L382.8,40.8 M473.8,128.9 L473.8,40.8 - M562.5,128.9 L562.5,40.8 '></path> -</g> -<g style="fill:none; color:blue; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter"> -</g> -<g style="fill:none; color:black; stroke:currentColor; stroke-width:1.00; stroke-linecap:butt; stroke-linejoin:miter"> - <path d='M93.4,40.8 L93.4,187.6 L582.7,187.6 L582.7,40.8 L93.4,40.8 Z '></path> -</g> -</svg> -
File [removed]: modes.dat Delta lines: +0 -8 =================================================================== --- examples/tov/modes.dat 2011-09-09 18:23:23 UTC (rev 130) +++ examples/tov/modes.dat 2011-09-13 04:29:02 UTC (rev 131) @@ -1,8 +0,0 @@ -1.44252691528028 F -3.95408396916149 H1 -5.91494894170517 H2 -7.77382493405710 H3 -9.58768293806276 H4 -11.3772129983097 H5 -13.1520241905666 H6 -14.9172321735655 H7