It is well known that python work wonders as a "glue language." Together with IPython magic functions we can execute code, push and pull variables, call functions, and even compile code from other languages in the IPython cell.
In previous posts I showed some examples with R, octave and matlab. So this time I will try a compiled language - Fortran. For that I will use the fantastic fortran_magic by MartÃn Gaitán.
Why Fortran? I will let the Simpsons explain why:
In [2]:
from IPython.display import YouTubeVideo
YouTubeVideo("0-hmzurm4SE")
Out[2]:
First we have to install the magic extension,
%install_ext https://raw.github.com/mgaitan/fortran_magic/master/fortranmagic.py
then load the extension:
In [3]:
%load_ext fortranmagic
MartÃn has a notebook with the documentation and examples. So I will go straight to my code and leave you to check his notebook.
The code I am compiling "interactively" is the MultiTaper Method by Mike Mann. Unless you are into climate/time-series analysis just go ahead and scroll down to the end of the code.
In [4]:
%%fortran
subroutine mtm_mann(spec_raw,spec_resh,spec_smoo,
$ spec_conf,recon_sig,time_serie,nscan,dt,npi,nwin,f1,f2,
$ inorm,ispec,iresh,ithresh,inoise,ilog,ismooth,
$ isignal,irecon,nsignals,iplotresh,iplotftest,iplotsmoo,
$ iplotraw,iplotconf,icon)
parameter (maxlen=50000,zero=0.0,maxsignal=100)
parameter (nlim=32768)
integer nscan,npi,nwin,inorm,ispec,iresh,ithresh
integer inoise,ilog,ismooth,isignal,irecon,nsignals,icon
integer iplotresh,iplotftest,iplotsmoo,iplotraw,iplotconf
real time_serie(nscan)
real spec_raw(nlim,3),spec_resh(nlim,3)
real spec_smoo(nlim,2),spec_conf(nlim,5)
real recon_sig(nlim,3)
real dt,f1,f2
cf2py intent(in) time_serie
cf2py intent(in) dt
cf2py intent(in) f1
cf2py intent(in) f2
cf2py intent(in) nscan
cf2py intent(hide) nscan
cf2py intent(in) npi
cf2py intent(in) nwin
cf2py intent(in) inorm,
cf2py intent(in) ispec
cf2py intent(in) iresh
cf2py intent(in) ithresh
cf2py intent(in) inoise
cf2py intent(in) ilog
cf2py intent(in) ismooth
cf2py intent(in) isignal
cf2py intent(in) irecon
cf2py intent(in) nsignals
cf2py intent(in) iplotresh
cf2py intent(in) iplotftest
cf2py intent(in) iplotsmoo
cf2py intent(in) iplotraw
cf2py intent(in) iplotconf
cf2py intent(in) icon
cf2py intent(out) spec_raw
cf2py intent(out) spec_smoo
cf2py intent(out) spec_conf
cf2py intent(out) spec_resh
cf2py intent(out) recon_sig
real a(maxlen),rcon(maxlen),signal(maxlen)
real fsignal(maxsignal),confsignal(maxsignal)
integer irec(maxsignal)
real anpi,fthresh,f0
integer nf0,nbnd
real fconf(6),conf(4)
real specraw0(nlim),specmed0(nlim),specresh0(nlim),
$ harmonic0(nlim),specback0(nlim),ftest(nlim)
real whiteraw0,whiterob0,rhoraw0,rhorob0,tauraw0,taurob0
do i=1,nscan
a(i) = time_serie(i)
end do
do i=nscan+1,maxlen
a(i)=zero
end do
demn=0
do i=1,nscan
demn = demn+a(i)
end do
demn = demn/float(nscan)
fray=1.0/(nscan*dt)
fny = 0.5/dt
anpi=float(npi)
fray=1.0/(nscan*dt)
fny = 0.5/dt
bndwdth = 2.0*anpi*fray
frange = f2-f1
if (ithresh.eq.0) fthresh=50.0
if (ithresh.eq.1) fthresh=90.0
if (ithresh.eq.2) fthresh=95.0
if (ithresh.eq.3) fthresh=99.0
if (ithresh.eq.4) fthresh=99.5
if (ithresh.eq.5) fthresh=99.9
fsmooth = frange/5.0
if (fsmooth.lt.bndwdth) fsmooth=bndwdth
do i=1,maxsignal
irec(i)=0
fsignal(i)=0.0
confsignal(i)=0.0
end do
call spec(a,nwin,npi,dt,nscan,f1,f2,
$ inoise,ismooth,ilog,fsmooth,
$ isignal,ispec,inorm,iresh,ithresh,
$ nf0,df,specraw0,specmed0,specresh0,harmonic0,
$ specback0,ftest,
$ conf,fconf,nbnd,
$ whiteraw0,whiterob0,rhoraw0,rhorob0,tauraw0,taurob0)
call findpeaks(flow,fhigh,isignal,
$ nf0,df,nbnd,specraw0,specmed0,specresh0,harmonic0,
$ specback0,ftest,
$ conf,fconf,
$ fsignal,confsignal,nsignals)
call display(nf0,f1,df,
$ spec_raw,spec_resh,spec_smoo,spec_conf,
$ specraw0,specmed0,specresh0,harmonic0,
$ specback0,ftest,conf,fconf,
$ whiteraw0,whiterob0,rhoraw0,rhorob0,
$ tauraw0,taurob0,
$ iplotresh,iplotftest,iplotsmoo,
$ iplotraw,iplotconf)
if (irecon.eq.1) then
do i=1,maxlen
signal(i)=0.0
rcon(i)=0.0
end do
do j=1,nsignals
if (irec(j).eq.1) then
f0 = fsignal(j)
call recon(a,nwin,npi,dt,nscan,bndwdth,f0,icon,rcon)
do i=1,nscan
signal(i)=signal(i)+rcon(i)
end do
endif
end do
do i=1,nscan
recon_sig(i,1)=float(i-1)*dt
recon_sig(i,2)=signal(i)
recon_sig(i,3)=a(i)-demn
end do
endif
end
subroutine spec(a0,nwin0,ip,dt0,nscan0,
$ flow,fhigh,
$ inoise0,ismooth,ilog0,fsmooth,
$ isignal,ispec,inorm,iresh,iper,
$ nf0,df,specraw0,specmed0,specresh0,harmonic0,specback0,
$ ftest,conf,fconf,nbnd,
$ whiteraw0,whiterob0,rhoraw0,rhorob0,tauraw0,taurob0)
integer maxwin,maxlen,maxdof,nlim,nmed
real big,small,tol
parameter (maxlen=50000)
parameter (maxwin=8,maxdof=2*maxwin,nlim=32768)
parameter (nmed=2000)
parameter (big=1.0e+36,small=1.0e-12,tol=1e-4)
character*2000 name21, name17
real dt0
real a0(maxlen)
integer ip,nscan0,nwin0,ilog0
integer inoise,inoise0
integer nbnd,nbnd0
integer idone(nlim)
real rho,rho0
real MSE
real*8 avar,el
real val(nmed)
real demean(nlim),specmed
real specraw0(nlim)
real specresh0(nlim),harmonic0(nlim),specback0(nlim),
$ specmed0(nlim)
real rednoise(nlim),rednoise0(nlim)
real ftest(nlim)
real base(nlim)
integer iharm(nlim)
real fconf(6)
real conf(4)
real adum,sumpeak
complex zed
real thresh
common/taperz/ta(16400,8),tai(16400,8),dcf(16400,8),
$ amu(16400,8)
common/staple/tad(32800),tad1(32800)
common/stap2/tas(32800,16)
common/misc/npad
common/npiprol/anpi
common/data/a(maxlen)
common/work/b(32800)
common/work2/fspec(16400)
common/BLK1/white,specmed(nlim),fnynew,ddf,f0
common/BLK2/ilog,nfreq
dimension dum(35)
real *8 b1(32800),b2(32800)
dimension reim(2),el(10)
equivalence (zed,reim),(iy,dum)
pi=3.14159265358979
radian=180./pi
tiny = 1e-6
nscan = nscan0
do i=1,maxlen
a(i)=a0(i)
end do
npi = ip
anpi=float(ip)
nwin=nwin0
dt = dt0
nscan = nscan0
ilog = ilog0
npts = nscan
npad=npts-1
ij=0
1000 npad=npad/2
ij=ij+1
if(npad.ne.0) go to 1000
npad=max(2**ij,1024)
fmin=flow
fmax=fhigh
fnynew = fhigh
ddf=1.0/(npad*dt)
df = ddf
eps=.1*ddf
n1=(fmin+eps)/ddf
n2=((fmax+eps)/ddf)
fmin=n1*ddf
fmax=n2*ddf
nf=n2-n1+1
nf0 = nf
nfreq=nf
f0 = fmin
nmin=n1*2+1
nmax=n2*2+2
fny = 0.5/dt
fray=1.0/(npts*dt)
bndwdth = 2.0*anpi*fray
halfwdth = bndwdth/2.0
nbnd = bndwdth/ddf
nbnd0 = fray/ddf
timefreq = 2.0*anpi
demn = 0.0
do i=1,nscan
demn=demn+a(i)
end do
demn = demn/float(nscan)
do i=1,nscan
demean(i)=a(i)-demn
end do
inoise = inoise0
if (inoise.gt.0) then
rho0 = 0.0
else
var = 0.0
do i=1,nscan
var = var + demean(i)**2
end do
var = var/float(nscan)
sd = sqrt(var)
c1 = 0.0
icount = 0
do i=2,nscan
icount = icount + 1
c1 = c1 + demean(i-1)*demean(i)
end do
c1 = c1/float(icount)
rho0 = c1/var
endif
name21='./data/chisquare.dat'
open (unit=21,file=name21,status='old')
do i=2,maxdof
idofs = 2*nwin
if (i.eq.idofs) then
read (21,*) idum,conf(1),conf(2),conf(3),conf(4)
else
read (21,*)
endif
end do
close (unit=21)
iftest = 1
jper = 3
if (nwin.eq.1) iftest = 0
if (iftest.eq.1) then
name17='./data/ftest.dat'
open (unit=17,file=name17,status='unknown')
do i=2,14,2
if (i.eq.2*nwin-2) then
read (17,*) idum,fconf(1),fconf(2),fconf(3),fconf(4),
$ fconf(5),fconf(6)
else
read (17,*)
endif
end do
close (unit=17)
thresh = fconf(iper)
endif
afact = 0.0
if ((iresh.eq.1).and.(isignal.ne.2)) then
afact = conf(1)
if (iper.eq.2) afact = conf(2)
if (iper.eq.3) afact = conf(3)
if (iper.gt.3) afact = conf(4)
endif
call taper2(npts,nwin,el)
if (inorm.eq.2) then ! amp spec (time limited)
anrm=1.0/dt
else
if(inorm.eq.1) then ! amp spect per unit time
anrm=float(npts)
else
anrm=1.
endif
endif
do iwin=1,nwin
do i=1,nscan
b1(i)=demean(i)*tas(i,iwin)
b2(i)=0.0
end do
j=nscan+1
do i=j,npad
b1(i)=0.0
b2(i)=0.0
end do
call fft2(b1,b2,npad)
sum=0.
j = 0
do i=1,npad/2
j = j+1
ta(j,iwin)=b1(i)/anrm
tai(j,iwin)=b2(i)/anrm
b(i) = b1(i)
b(i+1) = b2(i)
sumi=b1(i)*b1(i)+b2(i)*b2(i)
sum=sum+sumi
end do
end do
if (ispec.eq.0) then
call hires(ta,tai,el,nwin,nf-1,amu)
do i=1,nf-1
amu(i,1)=abs(amu(i,1))**2
end do
endif
if (ispec.eq.1) then
avar=0.d0
do i=1,nf-1
avar=avar+demean(i)*demean(i)
end do
if (inorm.eq.1) then ! amp spect per unit time
avar=avar/(npts*npts)
elseif(inorm.eq.2) then ! absolute amp spect
avar=avar*dt*dt
endif
call adwait(ta,tai,dcf,el,nwin,nf-1,amu,amu(1,2),avar)
do i=1,nf-1
amu(i,1)=abs(amu(i,1))
end do
endif
do i=1,nf-1
iharm(i)=0
specraw0(i)=amu(i,1)
specresh0(i)=amu(i,1)
harmonic0(i)=amu(i,1)
idone(i)=0
end do
if (iftest.eq.1) then
call regre(ta,tai,nf-1,nwin,amu)
do i=1,nf-1
ftest(i)=amu(i,3)
end do
endif
white0 = 0.0
do i=1,nf-1
white0 = white0+specraw0(i)
end do
white0 = white0/float(nf-1)
rho = rho0
white = white0
if (ismooth.eq.1) then
white = 0.0
nsmooth = fsmooth/ddf
do j=1,nf-1
if1 = j-nsmooth/2
if2 = j+nsmooth/2
if (if1.lt.1) then
if1=1
endif
if (if2.gt.nf-1) then
if2=nf-1
endif
nblk = if2-if1+1
do i=1,nblk
val(i)=specraw0(if1+i-1)
end do
kmid = (nblk+1)/2
do kk1=1,nblk
do kk2=kk1+1,nblk
if (val(kk2).lt.val(kk1)) then
adum = val(kk2)
val(kk2)=val(kk1)
val(kk1)=adum
endif
end do
end do
specmed(j)=val(kmid)
specmed0(j)=specmed(j)
white = white + specmed(j)
end do
white = white/float(nf-1)
if (inoise.ne.2) then
amin = big
do rho1=0.0,0.999,0.001
amiss = MSE(rho1)
if (amiss.lt.amin) then
rho=rho1
amin = amiss
endif
end do
endif
endif
tau = 1e+16
tau0 = 1e+16
if (rho.gt.0.0) tau = -dt/log(rho)
if (rho0.gt.0.0) tau0 = -dt/log(rho0)
do i=1,nf-1
ff = fmin+(i-1)*ddf
freq_norm = ff/fnynew
arg = freq_norm*pi
rednoise0(i) = white0*(1.0-rho0**2)/
$ (1.0-2.0*rho0*cos(arg)+rho0**2)
rednoise(i) = white*(1.0-rho**2)/
$ (1.0-2.0*rho*cos(arg)+rho**2)
if (isignal.eq.2) then
base(i)=0.0
else
if (inoise.lt.2) then
if (ismooth.eq.1) then
base(i) = rednoise(i)
else
base(i) = rednoise0(i)
endif
else
base(i) = specmed(i)*conf(1)
endif
specback0(i)=base(i)
endif
end do
if (iresh.eq.1) then
do i=nbnd,nf-1
if ((ftest(i).gt.thresh).
$ and.(specraw0(i).gt.afact*base(i))) iharm(i)=1
end do
do i=nbnd,nf-1
if (idone(i).ne.1) then
ipre0 = i-nbnd0/2-1
iaft0 = i+nbnd0/2+1
ipre = i-nbnd/2
iaft = i+nbnd/2
if (ipre.lt.nbnd+1) ipre=nbnd+1
if (iaft.gt.nf-2) iaft=nf-2
if (ipre0.lt.nbnd) ipre0=nbnd
if (iaft0.gt.nf-1) iaft0=nf-1
sumpeak = 0.0
if (iharm(i).eq.1) then
do j=ipre,iaft
specresh0(j)=0.5*(specraw0(ipre-1)+specraw0(iaft+1))
harmonic0(j)=specresh0(j)
end do
do j=ipre,iaft
idone(j)=1
sumpeak = sumpeak+specraw0(j)
end do
do j=ipre0,iaft0
harmonic0(j)=sumpeak/(iaft0-ipre0+1)
end do
endif
endif
end do
endif
if (isignal.eq.2) then
do i=1,nf-1
specback0(i)=specresh0(i)
end do
endif
whiteraw0=white0
whiterob0=white
rhoraw0=rho0
rhorob0=rho
tauraw0 = tau0
taurob0 = tau
8888 continue
return
end
real function MSE(rho)
parameter (nlim=32768)
COMMON /BLK1/white,specmed(nlim),fnynew,ddf,f0
COMMON /BLK2/ilog,nfreq
real rho
real pie,dff,freq_norm,ff,arg,small,rednoise,val1,val2
pie = 3.1415926
small = 1e-12
dff = 0.0
do j=1,nfreq
ff = f0+(j-1)*ddf
freq_norm = ff/fnynew
arg = freq_norm*pie
rednoise = white*(1.0-rho**2)/
$ (1.0-2.0*rho*cos(arg)+rho**2)
if (ilog.eq.0) then
dff = dff + (specmed(j)-rednoise)**2
else
val1 = abs(specmed(j))
val2 = abs(rednoise)
if (val1.lt.small) val1=small
if (val2.lt.small) val2=small
dff = dff +
$ (log(val1)-log(val2))**2
endif
end do
MSE = dff
return
end
subroutine hires(ta,tai,el,nwin,nf,ares)
real*8 el
dimension ta(16400,1),tai(16400,1),el(1),ares(1)
do j=1,nf
ares(j)=0.
end do
do i=1,nwin
a=1./(el(i)*nwin)
do j=1,nf
ares(j)=ares(j)+a*(ta(j,i)*ta(j,i)+tai(j,i)*tai(j,i))
end do
end do
do j=1,nf
ares(j)=sqrt(ares(j))
end do
return
end
subroutine adwait(ta,tai,dcf,el,nwin,nf,ares,degf,avar)
real*8 avar,spw,as,das,tol,el,a1,bias,scale,ax,fn,fx
dimension ta(16400,1),tai(16400,1),el(1),ares(1),degf(1)
dimension spw(10),bias(10),dcf(16400,1)
tol=3.d-4
jitter=0
scale=avar
do 200 i=1,nwin
200 bias(i)=(1.d0-el(i))
do 100 j=1,nf
do 150 i=1,nwin
150 spw(i)=(ta(j,i)*ta(j,i)+tai(j,i)*tai(j,i))/scale
as=(spw(1)+spw(2))/2.d0
do 300 k=1,20
fn=0.d0
fx=0.d0
do 350 i=1,nwin
a1=dsqrt(el(i))*as/(el(i)*as+bias(i))
a1=a1*a1
fn=fn+a1*spw(i)
fx=fx+a1
350 continue
ax=fn/fx
das=dabs(ax-as)
if(das/as.lt.tol) go to 400
300 as=ax
jitter=jitter+1
400 continue
ares(j)=as*scale
df=0.
do 450 i=1,nwin
dcf(j,i)=dsqrt(el(i))*as/(el(i)*as+bias(i))
450 df=df+dcf(j,i)*dcf(j,i)
100 degf(j)=df*2./(dcf(j,1)*dcf(j,1))
return
end
subroutine regre(sr,si,nf,nwin,amu)
real b
real *8 junk
dimension sr(16400,1),si(16400,1),amu(16400,1)
common/tapsum/b(10),junk(10)
sum=0.
do i=1,nwin
sum=sum+b(i)*b(i)
end do
do i=1,nf
amu(i,1)=0.
amu(i,2)=0.
do j=1,nwin
amu(i,1)=amu(i,1)+sr(i,j)*b(j)
amu(i,2)=amu(i,2)+si(i,j)*b(j)
end do
amu(i,1)=amu(i,1)/sum
amu(i,2)=amu(i,2)/sum
sum2=0.
do j=1,nwin
sumr=sr(i,j)-amu(i,1)*b(j)
sumi=si(i,j)-amu(i,2)*b(j)
sum2=sum2+sumr*sumr+sumi*sumi
end do
amu(i,3)=(nwin-1)*(amu(i,2)**2+amu(i,1)**2)*sum/sum2
end do
return
end
subroutine taper2(n,nwin,el)
real*8 el,a,z,pi,ww,cs,ai,an,eps,rlu,rlb
real*8 dfac,drat,gamma,bh,ell
real dump
common/npiprol/anpi
common/tapsum/tapsum(10),ell(10)
common/misc/npad
common/work/ip(32800)
common/taperz/z(65536),dump(393728)
common/stap2/ta(32800,16)
dimension a(32800,8),el(10)
data pi/3.14159265358979d0/
equivalence (a(1,1),ta(1,1))
an=dfloat(n)
ww=dble(anpi)/an
cs=dcos(2.d0*pi*ww)
do i=0,n-1
ai=dfloat(i)
a(i+1,1)=-cs*((an-1.d0)/2.d0-ai)**2
a(i+1,2)=-ai*(an-ai)/2.d0
a(i+1,3)=a(i+1,2)**2
end do
eps=1.e-13
m11=1
call tridib(n,eps,a(1,1),a(1,2),a(1,3),rlb,rlu,m11,nwin,el,ip,
x ierr,a(1,4),a(1,5))
call tinvit(n,n,a(1,1),a(1,2),a(1,3),nwin,el,ip,z,ierr,
x a(1,4),a(1,5),a(1,6),a(1,7),a(1,8))
dfac=an*pi*ww
drat=8.d0*dfac
dfac=4.d0*dsqrt(pi*dfac)*dexp(-2.d0*dfac)
do k=1,nwin
el(k)=1.d0-dfac
dfac=dfac*drat/k ! is this correct formula? yes,but fails as k -> 2n
end do
gamma=dlog(8.d0*an*dsin(2.d0*pi*ww))+0.5772156649d0
do k=1,nwin
bh=-2.d0*pi*(an*ww-dfloat(k-1)/2.d0-.25d0)/gamma
ell(k)=1.d0/(1.d0+dexp(pi*bh))
end do
do i=1,nwin
el(i)=dmax1(ell(i),el(i))
end do
do k=1,nwin
kk=(k-1)*n
tapsum(k)=0.
tapsq=0.
do i=1,n
aa=z(kk+i)
ta(i,k)=aa
tapsum(k)=tapsum(k)+aa
tapsq=tapsq+aa*aa
end do
aa=sqrt(tapsq/n)
tapsum(k)=tapsum(k)/aa
do i=1,n
ta(i,k)=ta(i,k)/aa
end do
end do
return
end
subroutine findpeaks(flow,fhigh,isignal,
$ nf0,df,nbnd,specraw0,specmed0,specresh0,harmonic0,
$ specback0,ftest,
$ conf,fconf,
$ fsignal,confsignal,nsignals)
parameter (maxsignal=40,nlarge=300)
parameter (nlim=32768)
real fsignal(maxsignal),confsignal(maxsignal)
real fsig(nlarge),confsig(nlarge),sigrat(nlarge)
integer null(nlim)
integer nbnd,isignal
real fconf(6),conf(4),sig(3)
real ratio,ratmax,ratio0,dummy
real specraw0(nlim),specmed0(nlim),specresh0(nlim),
$ harmonic0(nlim),specback0(nlim),ftest(nlim)
nsignals = 0
do i=1,nf0
null(i)=0
end do
sig(1)=99.0
sig(2)=95.0
sig(3)=90.0
do k=1,3
if (isignal.lt.2) then
thresh = conf(5-k)
else
thresh = fconf(5-k)
endif
do i=1,nf0-1
ff = flow+float(i-1)*df
if (isignal.lt.2) then
ratio0 = specraw0(i)/specback0(i)
else
ratio0 = ftest(i)
endif
if ((ratio0.gt.thresh).and.(null(i).eq.0)) then
nsignals=nsignals+1
if (nsignals.gt.nlarge) then
nsignals=nsignals-1
goto 888
endif
ratmax = ratio0
imax = i
fmax = ff
do j=max(1,i-nbnd),min(i+nbnd,nf0-1)
f0 = flow+float(j-1)*df
if (isignal.lt.2) then
ratio = specraw0(j)/specback0(j)
else
ratio = ftest(j)
endif
if (ratio.gt.ratmax) then
imax = j
ratmax=ratio
fmax = f0
endif
end do
if (null(imax).eq.1) then
nsignals = nsignals-1
goto 444
endif
fsig(nsignals)=fmax
confsig(nsignals)=sig(k)
sigrat(nsignals)=ratmax
do j=max(1,imax-nbnd),min(imax+nbnd,nf0-1)
null(j)=1
end do
endif
444 continue
end do
end do
888 continue
do i=1,nsignals
do j=i+1,nsignals
if (sigrat(j).gt.sigrat(i)) then
dummy = sigrat(j)
sigrat(j)=sigrat(i)
sigrat(i) = dummy
dummy = fsig(j)
fsig(j)=fsig(i)
fsig(i) = dummy
dummy = confsig(j)
confsig(j)=confsig(i)
confsig(i) = dummy
endif
end do
end do
nsignals=min(maxsignal,nsignals)
do i=1,nsignals
do j=i+1,nsignals
if (fsig(j).lt.fsig(i)) then
dummy = fsig(i)
fsig(i)=fsig(j)
fsig(j)=dummy
dummy = confsig(i)
confsig(i)=confsig(j)
confsig(j)=dummy
endif
end do
end do
do i=1,nsignals
fsignal(i)=fsig(i)
confsignal(i)=confsig(i)
end do
return
end
subroutine display(nf0,f1,df,
$ spec_raw,spec_resh,spec_smoo,spec_conf,
$ specraw0,specmed0,specresh0,harmonic0,
$ specback0,ftest,conf,fconf,
$ whiteraw0,whiterob0,rhoraw0,rhorob0,
$ tauraw0,taurob0,
$ iplotresh,iplotftest,iplotsmoo,
$ iplotraw,iplotconf)
parameter (nlim=32768)
integer iplotraw,iplotresh,iplotftest,iplotsmoo,
$ iplotconf
real fconf(6),conf(4)
real spec_raw(nlim,3),spec_resh(nlim,3)
real spec_smoo(nlim,2),spec_conf(nlim,5)
real specraw0(nlim),specmed0(nlim),specresh0(nlim),
$ harmonic0(nlim),specback0(nlim),ftest(nlim)
real whiteraw0,whiterob0,rhoraw0,rhorob0,tauraw0,taurob0
real f1,df
integer nf0
do i=1,nf0-1
ff = f1+float(i-1)*df
if (iplotraw.eq.1) then
spec_raw(i,1) = ff
spec_raw(i,2) = specraw0(i)
spec_raw(i,3) = ftest(i)
endif
if (iplotresh.eq.1) then
spec_resh(i,1) = ff
spec_resh(i,2) = specresh0(i)
spec_resh(i,3) = harmonic0(i)
endif
if (iplotsmoo.eq.1) then
spec_smoo(i,1) = ff
spec_smoo(i,2) = specmed0(i)
endif
if (iplotconf.eq.1) then
spec_conf(i,1) = ff
spec_conf(i,2) = specback0(i)*conf(1)
spec_conf(i,3) = specback0(i)*conf(2)
spec_conf(i,4) = specback0(i)*conf(3)
spec_conf(i,5) = specback0(i)*conf(4)
endif
end do
8888 continue
652 format (f7.4,e32.8)
653 format (f7.4,2e30.8)
654 format (f7.4,1x,4e17.3)
655 format (f7.4,1x,5e14.2)
return
end
subroutine recon(a0,nwin0,npi,dt0,nscan,bndwdth,ff0,icon,rcon)
parameter (maxlen=50000,maxsignal=20)
real a0(maxlen),dt0
integer nwin0,nscan,npi
real*8 ell,env0,cs,sn,ex
real bndwdth,f0,ff0
character*4 chead(158)
complex zed
common/npiprol/anpi
common/header/iah(158)
common/data/a(maxlen)
common/nnaamm/iwflag
common/work/b(8194)
common/pltopt/npad,nst,mpts,nmin,nmax,t(3),fmin,fmax,nf,
x ddt,nwin,inc,tt(3)
common/taperz/ta(4100,16),tai(4100,16)
common/staple/tad(4000)
common/stap2/tas(8192,16),ex(500)
common/stap4/cs(8200),sn(8200)
common/envel1/ell(20)
dimension dum(35),ahead(158),env0(2,1000)
dimension reim(2),time(2)
real *8 b1(8194),b2(8194)
real env(3,500,2),aa(8200,3),acent(8200),rcon(maxlen)
equivalence (zed,reim),(iy,dum),(iah,ahead),(iah,chead),
x (tad,env0)
iwflag=0
pi=3.14159265358979
do i=1,maxlen
a(i)=a0(i)
end do
time(1)=0.
dt = dt0
time(2)=dt/1.
fzero=0.
fny=0.5/dt
npts=nscan
nwin = nwin0
anpi = float(npi)
f1 = 0.0
f2 = fny
npad=2
do while (npad.lt.nscan)
npad=npad*2
end do
npad=8192
ddf=1./(npad*dt)
eps=.1*ddf
nn1=(f1+eps)/ddf
fmin=nn1*ddf
nn2=(f2+eps)/ddf
fmax=nn2*ddf
nf=nn2-nn1+1
nmin=nn1*2+1
nmax=nn2*2+2
t(1)=fmin
t(2)=ddf
t(3)=ddf
ftfrq=ddf
102 format(i6,g12.4,i6)
n1=1
n2=nscan
2660 call taper1(npts,nwin,ell)
demn=0.
do i=n1,n2
demn=demn+a(i)
end do
demn=demn/(n2-n1+1)
do i=n1,n2
acent(i)=a(i)-demn
end do
anrm=npts
icont=0
do iwin=1,nwin
do i=n1,n2
b1(i-n1+1)=acent(i)*tas(i+1-n1,iwin)
b2(i-n1+1)=0.0
end do
j=npts+1
do i=j,npad
b1(i)=0.
b2(i)=0.
end do
call fft2(b1,b2,npad)
j = 0
do i=1,npad/2
j = j+1
b(i)=b1(i)
b(i+1)=b2(i)
ta(i,iwin)=b1(i)
tai(i,iwin)=b2(i)
end do
end do
1660 inc=(npts-1)/500+1
j=0
do i=1,npts,inc
j=j+1
do k=1,nwin
tas(j,k)=tas(i,k)
end do
end do
mpts=j
ddt=dt*inc
f0=ff0
if (f0.lt.bndwdth) f0=0.0
anorm = 1.0
if (f0.lt.eps) anorm=0.5
jmax = 1
if (icon.eq.0) jmax=3
do i=1,jmax
itype = icon-1
if (icon.eq.0) itype=i-1
call envel(f0,itype)
do j=1,mpts
env(i,j,1)=env0(1,j)
env(i,j,2)=env0(2,j)
end do
end do
call cossin(npts,cs,sn,1.d0,0.d0,dble(f0),dble(dt))
finc=float(inc)
do jj=1,jmax
do i=1,mpts
env0(1,i)=env(jj,i,1)
env0(2,i)=env(jj,i,2)
end do
do i=1,mpts-1
ii=(i-1)*inc
do j=1,inc
a1=(env0(1,i)*(inc+1-j)+env0(1,i+1)*(j-1))/finc
a2=(env0(2,i)*(inc+1-j)+env0(2,i+1)*(j-1))/finc
a(n1-1+nscan+ii+j)=a1*cs(ii+j)+a2*sn(ii+j)
aa(n1-1+ii+j,jj)=anorm*a(n1-1+nscan+ii+j)
end do
end do
ii=(mpts-1)*inc
e1=2.*env0(1,mpts)-env0(1,mpts-1)
e2=2.*env0(2,mpts)-env0(2,mpts-1)
do j=1,inc
a1=(env0(1,mpts)*(inc+1-j)+e1*(j-1))/finc
a2=(env0(2,mpts)*(inc+1-j)+e2*(j-1))/finc
a(n1-1+nscan+ii+j)=a1*cs(ii+j)+a2*sn(ii+j)
aa(n1-1+ii+j,jj)=anorm*a(n1-1+nscan+ii+j)
end do
end do
w1keep = 1.0
w2keep = 0.0
w3keep = 0.0
abest = 1.0e+36
step = 0.02
if (icon.eq.0) then
do w1=0,1,step
do w2=0,1,step
w3 = 1.0-w1-w2
if ((w3.ge.0.0).and.(w3.le.1.0)) then
amiss = 0.0
do i=1,nscan
signal = w1*aa(i,1)+w2*aa(i,2)+w3*aa(i,3)
amiss = amiss + (signal-acent(i))**2
end do
if (amiss.lt.abest) then
abest=amiss
w1keep = w1
w2keep = w2
w3keep = w3
endif
endif
end do
end do
do i=1,nscan
rcon(i)=w1keep*aa(i,1)+w2keep*aa(i,2)+w3keep*aa(i,3)
end do
else
do i=1,nscan
rcon(i)=aa(i,1)
end do
do i=nscan+1,maxlen
rcon(i)=0.0
end do
endif
return
end
subroutine taper1(n,nwin,el)
real*8 el,a,z,pi,ww,cs,ai,an,eps,rlu,rlb,ex
real*8 dfac,drat,gamma,bh,ell
common/nnaamm/iwflag
common/npiprol/anpi
common/tapsum/tapsum(20),ell(20)
common/work/ip(8194)
common/taperz/z(65600) ! we use this common block for scratch space
common/stap2/ta(8192,16),ex(500)
dimension a(8192,8),el(10)
data iwflag/0/,pi/3.14159265358979d0/
equivalence (a(1,1),ta(1,1))
an=dfloat(n)
ww=dble(anpi)/an
cs=dcos(2.d0*pi*ww)
do i=0,n-1
ai=dfloat(i)
a(i+1,1)=-cs*((an-1.d0)/2.d0-ai)**2
a(i+1,2)=-ai*(an-ai)/2.d0
a(i+1,3)=a(i+1,2)**2 ! required by eispack routine
end do
eps=1.e-13
m11=1
call tridib(n,eps,a(1,1),a(1,2),a(1,3),rlb,rlu,m11,nwin,el,ip,
x ierr,a(1,4),a(1,5))
call tinvit(n,n,a(1,1),a(1,2),a(1,3),nwin,el,ip,z,ierr,
x a(1,4),a(1,5),a(1,6),a(1,7),a(1,8))
dfac=an*pi*ww
drat=8.d0*dfac
dfac=4.d0*dsqrt(pi*dfac)*dexp(-2.d0*dfac)
do k=1,nwin
el(k)=1.d0-dfac
dfac=dfac*drat/k ! note that this doesn't hold as k -> 2n
end do
gamma=dlog(8.d0*an*dsin(2.d0*pi*ww))+0.5772156649d0
do k=1,nwin
bh=-2.d0*pi*(an*ww-dfloat(k-1)/2.d0-.25d0)/gamma
ell(k)=1.d0/(1.d0+dexp(pi*bh))
end do
do i=1,nwin
el(i)=dmax1(ell(i),el(i))
end do
do k=1,nwin
kk=(k-1)*n
tapsum(k)=0.
tapsq=0.
do i=1,n
aa=z(kk+i)
ta(i,k)=aa
tapsum(k)=tapsum(k)+aa
tapsq=tapsq+aa*aa
end do
aa=sqrt(tapsq/n)
tapsum(k)=tapsum(k)/aa
do i=1,n
ta(i,k)=ta(i,k)/aa
end do
end do
101 format(80a)
iwflag=1
return
end
subroutine envel(f0,inorm)
real*8 cs,sn,env0,amp,ell,dex,ex
real qc
complex*16 g,d,env,dum,za,amp0,amp1,qrsave
common/pltopt/npad,nst,npts,nmin,nmax,t(3),fmin,fmax,nf,
x dt,nwin,inc,tt(3)
common/envel1/ell(20)
common/taperz/ta(4100,16),tai(4100,16)
common/staple/env(1000)
common/stap2/tas(8192,16),ex(500)
common/stap3/g(150000)
common/stap4/cs(8200),sn(8200)
common/stap5/d(100),za(100),qrsave(100),dum2(2000),dum(8200)
dimension env0(2,1000),amp(2)
equivalence (env0,env),(amp0,amp)
tt(1)=0.
tt(2)=dt
pi=3.14159265358979
rad=180./pi
eps=1.e-4
jsmoo = inorm
iboost=0
fmin=t(1)
fmax=t(1)+(nf-1)*t(2)
ff0 = f0
qc=0.
dex=dble(qc*pi/npts)
call boost(iboost,npts,ex,dex)
call fillup(f0)
call premult(iboost,jsmoo)
call zqrdc(g,npts,npts,nwin,qrsave,ip,dumm,0)
call zbacktr(nwin,npts,g,d,dum,2)
amp0=dcmplx(0.d0,0.d0)
do k=1,nwin
amp0=amp0+conjg(za(k))*dum(k) ! note: conjugated
end do
call zbacktr(nwin,npts,g,za,dum,2)
amp1=dcmplx(0.d0,0.d0)
do k=1,nwin
amp1=amp1+conjg(za(k))*dum(k) ! note: conjugated
end do
amp0=amp0/amp1
sum1=0.
do k=1,nwin
sum1=sum1+(cdabs(d(k)))**2
end do
do k=1,nwin
d(k)=d(k)-za(k)*amp0
end do
sum2=0.
do k=1,nwin
sum2=sum2+(cdabs(d(k)))**2
end do
rat=sum2/sum1
do i=1,npts
env(i)=dcmplx(0.,0.)
end do
call zbacktr(nwin,npts,g,d,env,1)
call zqrsl(g,npts,npts,nwin,qrsave,env,env,
x dum,dum,dum,dum,10000,info)
call postmult(jsmoo,npts,env,amp0,cs)
npts2=npts*2
return
end
subroutine zbacktr(n,npts,g,d,dum,ick)
complex*16 g(npts,1),dum(1),d(1)
do k=1,n
dum(k)=d(k)
end do
do i=1,n
i1=i-1
if(i.gt.1) then
do j=1,i1
dum(i)=dum(i)-conjg(g(j,i))*dum(j)
end do
endif
dum(i)=dum(i)/conjg(g(i,i))
end do
if(ick.eq.1) return
do i=n,1,-1
ip1=i+1
if(i.lt.n) then
do j=ip1,n
dum(i)=dum(i)-g(i,j)*dum(j)
end do
endif
dum(i)=dum(i)/g(i,i)
end do
return
end
subroutine boost(iboost,npts,ex,dex)
implicit real*8 (a-h,o-z)
dimension ex(1)
if(iboost.eq.1) then
dex=dexp(-dex)
ex(1)=1.d0
do i=2,npts
ex(i)=ex(i-1)*dex
end do
else ! dont boost the envelope
do i=1,npts
ex(i)=1.d0
end do
endif
return
end
subroutine fillup(f0)
real*8 cs,sn,ex
complex*16 g,d,junk1,junk3
real *4 junk2
common/pltopt/npad,nst,npts,nmin,nmax,t(3),fmin,fmax,nf,
x dt,nwin,inc,tt(3)
common/taperz/ta(4100,16),tai(4100,16)
common/stap2/tas(8192,16),ex(500)
common/stap3/g(150000)
common/stap4/cs(8200),sn(8200)
common/stap5/d(100),junk1(200),junk2(2000),junk3(8200)
n1=(f0-fmin)/t(2)+1
f1=(fmin+(n1-1)*t(2))
df1=f0-f1
if(df1*2.0.gt.t(2)) then
n1=n1+1
f1=f1+t(2)
df1=f1-f0
else
df1=-df1 ! since formula uses f1-f0
endif
call cossin(npts,cs,sn,1.d0,0.d0,dble(df1),dble(dt))
do k=1,nwin
d(k)=dcmplx(ta(n1,k),tai(n1,k))*2.d0/inc
end do
do k=1,nwin
ii=(k-1)*npts
jj=(k-1)*npts
j=0
do i=1,npts
j=j+1
g(jj+j)=ex(i)*dcmplx(tas(i,k)*cs(j),-tas(i,k)*sn(j))
end do
end do
return
end
subroutine premult(iboost,jsmoo)
real*8 cs,sn,cs0
complex*16 g,d,za,junk1,junk3
real *4 junk2
common/pltopt/npad,nst,npts,nmin,nmax,t(3),fmin,fmax,nf,
x dt,nwin,inc,tt(3)
common/stap3/g(150000)
common/stap4/cs(8200),sn(8200)
common/stap5/d(100),za(100),junk1(100),
x junk2(2000),junk3(8200)
pi=3.14159265358979
eps=1.e-4
if(iboost.eq.0) then
cs0=qc*pi/(npts-1)
cs0=dexp(-cs0)
cs(1)=1.d0
do i=2,npts
cs(i)=cs(i-1)*cs0
end do
else ! if we have boosted the kernels, the fixed element is constant
do i=1,npts
cs(i)=1.d0
end do
endif
do k=1,nwin
za(k)=dcmplx(0.,0.)
kk=(k-1)*npts
do i=1,npts
za(k)=za(k)+cs(i)*g(kk+i)
end do
za(k)=conjg(za(k)) ! must conjugate to adhere to conventions
end do
if(jsmoo.gt.0) then
do ksmoo=1,jsmoo
do k=1,nwin
ii=(k-1)*npts
nptsm=npts-1
do i=nptsm,1,-1
g(ii+i)=g(ii+i)+g(ii+i+1)
end do
g(ii+1)=g(ii+1)/eps
end do
end do
endif
return
end
subroutine postmult(jsmoo,npts,env,amp0,cs)
real*8 cs(1)
complex*16 env(1),amp0
eps=1.e-4
if(jsmoo.gt.0) then
do ksmoo=1,jsmoo
env(1)=env(1)/eps
do i=2,npts
env(i)=env(i)+env(i-1)
end do
end do
endif
do i=1,npts
env(i)=env(i)+amp0*cs(i)
end do
return
end
subroutine gfill(npts,nwin,ib,jb,df,ddt,g,ex,tas)
implicit real*8 (a-h,o-z)
real*4 tas,tt(3)
real*8 junk1
complex*16 g
common/stap4/cs(500),sn(500),junk1(15400)
dimension g(npts,2,nwin,2),ex(1),tas(8192,1)
call cossin(npts,cs,sn,1.d0,0.d0,df,ddt)
tt(1)=0.
tt(2)=ddt
do k=1,nwin
do i=1,npts
g(i,ib,k,jb)=ex(i)*dcmplx(tas(i,k)*cs(i),-tas(i,k)*sn(i))
end do
end do
return
end
SUBROUTINE COSSIN(N,C,S,C0,S0,F,DT)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION C(1),S(1)
DATA PI/3.14159265358979D0/
C(1)=C0
S(1)=S0
CS=DCOS(2.D0*PI*F*DT)
SN=DSIN(2.D0*PI*F*DT)
DO 100 I=2,N
C(I)=C(I-1)*CS-S(I-1)*SN
100 S(I)=C(I-1)*SN+S(I-1)*CS
RETURN
END
SUBROUTINE FFT2(ar,ai,N)
implicit real*8 (a-h,o-z)
implicit integer*4 (i-n)
dimension ar(1),ai(1)
mex=dlog(dble(float(n)))/.693147d0
nv2=N/2
nm1=N-1
j=1
do 7 i=1,nm1
if(i .ge. j) go to 5
tr=ar(j)
ti=ai(j)
ar(j)=ar(i)
ai(j)=ai(i)
ar(i)=tr
ai(i)=ti
5 k=nv2
6 if(k .ge. j) go to 7
j=j-k
k=k/2
go to 6
7 j=j+k
pi=3.14159265358979d0
do 20 l=1,mex
le=2**l
le1=le/2
ur=1.0
ui=0.
wr=dcos(pi/le1 )
wi=dsin (pi/le1)
do 20 j=1,le1
do 10 i=j,N,le
ip=i+le1
tr=ar(ip)*ur - ai(ip)*ui
ti=ai(ip)*ur + ar(ip)*ui
ar(ip)=ar(i)-tr
ai(ip)=ai(i) - ti
ar(i)=ar(i)+tr
ai(i)=ai(i)+ti
10 continue
utemp=ur
ur=ur*wr - ui*wi
ui=ui*wr + utemp*wi
20 continue
return
end
SUBROUTINE TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,M,W,IND,IERR,RV4,RV5)
INTEGER I,J,K,L,M,N,P,Q,R,S,II,M1,M2,M11,M22,TAG,IERR,ISTURM
DOUBLE PRECISION D(N),E(N),E2(N),W(M),RV4(N),RV5(N)
DOUBLE PRECISION U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,TST1,TST2,EPSLON
INTEGER IND(M)
IERR = 0
TAG = 0
XU = D(1)
X0 = D(1)
U = 0.0D0
DO 40 I = 1, N
X1 = U
U = 0.0D0
IF (I .NE. N) U = DABS(E(I+1))
XU = DMIN1(D(I)-(X1+U),XU)
X0 = DMAX1(D(I)+(X1+U),X0)
IF (I .EQ. 1) GO TO 20
TST1 = DABS(D(I)) + DABS(D(I-1))
TST2 = TST1 + DABS(E(I))
IF (TST2 .GT. TST1) GO TO 40
20 E2(I) = 0.0D0
40 CONTINUE
X1 = N
X1 = X1 * EPSLON(DMAX1(DABS(XU),DABS(X0)))
XU = XU - X1
T1 = XU
X0 = X0 + X1
T2 = X0
P = 1
Q = N
M1 = M11 - 1
IF (M1 .EQ. 0) GO TO 75
ISTURM = 1
50 V = X1
X1 = XU + (X0 - XU) * 0.5D0
IF (X1 .EQ. V) GO TO 980
GO TO 320
60 IF (S - M1) 65, 73, 70
65 XU = X1
GO TO 50
70 X0 = X1
GO TO 50
73 XU = X1
T1 = X1
75 M22 = M1 + M
IF (M22 .EQ. N) GO TO 90
X0 = T2
ISTURM = 2
GO TO 50
80 IF (S - M22) 65, 85, 70
85 T2 = X1
90 Q = 0
R = 0
100 IF (R .EQ. M) GO TO 1001
TAG = TAG + 1
P = Q + 1
XU = D(P)
X0 = D(P)
U = 0.0D0
DO 120 Q = P, N
X1 = U
U = 0.0D0
V = 0.0D0
IF (Q .EQ. N) GO TO 110
U = DABS(E(Q+1))
V = E2(Q+1)
110 XU = DMIN1(D(Q)-(X1+U),XU)
X0 = DMAX1(D(Q)+(X1+U),X0)
IF (V .EQ. 0.0D0) GO TO 140
120 CONTINUE
140 X1 = EPSLON(DMAX1(DABS(XU),DABS(X0)))
IF (EPS1 .LE. 0.0D0) EPS1 = -X1
IF (P .NE. Q) GO TO 180
IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940
M1 = P
M2 = P
RV5(P) = D(P)
GO TO 900
180 X1 = X1 * (Q - P + 1)
LB = DMAX1(T1,XU-X1)
UB = DMIN1(T2,X0+X1)
X1 = LB
ISTURM = 3
GO TO 320
200 M1 = S + 1
X1 = UB
ISTURM = 4
GO TO 320
220 M2 = S
IF (M1 .GT. M2) GO TO 940
X0 = UB
ISTURM = 5
DO 240 I = M1, M2
RV5(I) = UB
RV4(I) = LB
240 CONTINUE
K = M2
250 XU = LB
DO 260 II = M1, K
I = M1 + K - II
IF (XU .GE. RV4(I)) GO TO 260
XU = RV4(I)
GO TO 280
260 CONTINUE
280 IF (X0 .GT. RV5(K)) X0 = RV5(K)
300 X1 = (XU + X0) * 0.5D0
IF ((X0 - XU) .LE. DABS(EPS1)) GO TO 420
TST1 = 2.0D0 * (DABS(XU) + DABS(X0))
TST2 = TST1 + (X0 - XU)
IF (TST2 .EQ. TST1) GO TO 420
320 S = P - 1
U = 1.0D0
DO 340 I = P, Q
IF (U .NE. 0.0D0) GO TO 325
V = DABS(E(I)) / EPSLON(1.0D0)
IF (E2(I) .EQ. 0.0D0) V = 0.0D0
GO TO 330
325 V = E2(I) / U
330 U = D(I) - X1 - V
IF (U .LT. 0.0D0) S = S + 1
340 CONTINUE
GO TO (60,80,200,220,360), ISTURM
360 IF (S .GE. K) GO TO 400
XU = X1
IF (S .GE. M1) GO TO 380
RV4(M1) = X1
GO TO 300
380 RV4(S+1) = X1
IF (RV5(S) .GT. X1) RV5(S) = X1
GO TO 300
400 X0 = X1
GO TO 300
420 RV5(K) = X1
K = K - 1
IF (K .GE. M1) GO TO 250
900 S = R
R = R + M2 - M1 + 1
J = 1
K = M1
DO 920 L = 1, R
IF (J .GT. S) GO TO 910
IF (K .GT. M2) GO TO 940
IF (RV5(K) .GE. W(L)) GO TO 915
DO 905 II = J, S
I = L + S - II
W(I+1) = W(I)
IND(I+1) = IND(I)
905 CONTINUE
910 W(L) = RV5(K)
IND(L) = TAG
K = K + 1
GO TO 920
915 J = J + 1
920 CONTINUE
940 IF (Q .LT. N) GO TO 100
GO TO 1001
980 IERR = 3 * N + ISTURM
1001 LB = T1
UB = T2
RETURN
END
SUBROUTINE TINVIT(NM,N,D,E,E2,M,W,IND,Z,
X IERR,RV1,RV2,RV3,RV4,RV6)
INTEGER I,J,M,N,P,Q,R,S,II,IP,JJ,NM,ITS,TAG,IERR,GROUP
DOUBLE PRECISION D(N),E(N),E2(N),W(M),Z(NM,M),
X RV1(N),RV2(N),RV3(N),RV4(N),RV6(N)
DOUBLE PRECISION U,V,UK,XU,X0,X1,EPS2,EPS3,EPS4,NORM,ORDER,EPSLON,
X PYTHAG
INTEGER IND(M)
IERR = 0
IF (M .EQ. 0) GO TO 1001
TAG = 0
ORDER = 1.0D0 - E2(1)
Q = 0
100 P = Q + 1
DO 120 Q = P, N
IF (Q .EQ. N) GO TO 140
IF (E2(Q+1) .EQ. 0.0D0) GO TO 140
120 CONTINUE
140 TAG = TAG + 1
S = 0
DO 920 R = 1, M
IF (IND(R) .NE. TAG) GO TO 920
ITS = 1
X1 = W(R)
IF (S .NE. 0) GO TO 510
XU = 1.0D0
IF (P .NE. Q) GO TO 490
RV6(P) = 1.0D0
GO TO 870
490 NORM = DABS(D(P))
IP = P + 1
DO 500 I = IP, Q
500 NORM = DMAX1(NORM, DABS(D(I))+DABS(E(I)))
EPS2 = 1.0D-3 * NORM
EPS3 = EPSLON(NORM)
UK = Q - P + 1
EPS4 = UK * EPS3
UK = EPS4 / DSQRT(UK)
S = P
505 GROUP = 0
GO TO 520
510 IF (DABS(X1-X0) .GE. EPS2) GO TO 505
GROUP = GROUP + 1
IF (ORDER * (X1 - X0) .LE. 0.0D0) X1 = X0 + ORDER * EPS3
520 V = 0.0D0
DO 580 I = P, Q
RV6(I) = UK
IF (I .EQ. P) GO TO 560
IF (DABS(E(I)) .LT. DABS(U)) GO TO 540
XU = U / E(I)
RV4(I) = XU
RV1(I-1) = E(I)
RV2(I-1) = D(I) - X1
RV3(I-1) = 0.0D0
IF (I .NE. Q) RV3(I-1) = E(I+1)
U = V - XU * RV2(I-1)
V = -XU * RV3(I-1)
GO TO 580
540 XU = E(I) / U
RV4(I) = XU
RV1(I-1) = U
RV2(I-1) = V
RV3(I-1) = 0.0D0
560 U = D(I) - X1 - XU * V
IF (I .NE. Q) V = E(I+1)
580 CONTINUE
IF (U .EQ. 0.0D0) U = EPS3
RV1(Q) = U
RV2(Q) = 0.0D0
RV3(Q) = 0.0D0
600 DO 620 II = P, Q
I = P + Q - II
RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
V = U
U = RV6(I)
620 CONTINUE
IF (GROUP .EQ. 0) GO TO 700
J = R
DO 680 JJ = 1, GROUP
630 J = J - 1
IF (IND(J) .NE. TAG) GO TO 630
XU = 0.0D0
DO 640 I = P, Q
640 XU = XU + RV6(I) * Z(I,J)
DO 660 I = P, Q
660 RV6(I) = RV6(I) - XU * Z(I,J)
680 CONTINUE
700 NORM = 0.0D0
DO 720 I = P, Q
720 NORM = NORM + DABS(RV6(I))
IF (NORM .GE. 1.0D0) GO TO 840
IF (ITS .EQ. 5) GO TO 830
IF (NORM .NE. 0.0D0) GO TO 740
RV6(S) = EPS4
S = S + 1
IF (S .GT. Q) S = P
GO TO 780
740 XU = EPS4 / NORM
DO 760 I = P, Q
760 RV6(I) = RV6(I) * XU
780 DO 820 I = IP, Q
U = RV6(I)
IF (RV1(I-1) .NE. E(I)) GO TO 800
U = RV6(I-1)
RV6(I-1) = RV6(I)
800 RV6(I) = U - RV4(I) * RV6(I-1)
820 CONTINUE
ITS = ITS + 1
GO TO 600
830 IERR = -R
XU = 0.0D0
GO TO 870
840 U = 0.0D0
DO 860 I = P, Q
860 U = PYTHAG(U,RV6(I))
XU = 1.0D0 / U
870 DO 880 I = 1, N
880 Z(I,R) = 0.0D0
DO 900 I = P, Q
900 Z(I,R) = RV6(I) * XU
X0 = X1
920 CONTINUE
IF (Q .LT. N) GO TO 100
1001 RETURN
END
DOUBLE PRECISION FUNCTION EPSLON (X)
DOUBLE PRECISION X
DOUBLE PRECISION A,B,C,EPS
A = 4.0D0/3.0D0
10 B = A - 1.0D0
C = B + B + B
EPS = DABS(C-1.0D0)
IF (EPS .EQ. 0.0D0) GO TO 10
EPSLON = EPS*DABS(X)
RETURN
END
DOUBLE PRECISION FUNCTION PYTHAG(A,B)
DOUBLE PRECISION A,B
DOUBLE PRECISION P,R,S,T,U
P = DMAX1(DABS(A),DABS(B))
IF (P .EQ. 0.0D0) GO TO 20
R = (DMIN1(DABS(A),DABS(B))/P)**2
10 CONTINUE
T = 4.0D0 + R
IF (T .EQ. 4.0D0) GO TO 20
S = R/T
U = 1.0D0 + 2.0D0*S
P = U*P
R = (S/U)**2 * R
GO TO 10
20 PYTHAG = P
RETURN
END
SUBROUTINE ZQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,JOB)
INTEGER LDX,N,P,JOB
INTEGER JPVT(1)
COMPLEX*16 X(LDX,1),QRAUX(1),WORK(1)
INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU
DOUBLE PRECISION MAXNRM,DZNRM2,TT
COMPLEX*16 ZDOTC,NRMXL,T
LOGICAL NEGJ,SWAPJ
COMPLEX*16 CSIGN,ZDUM,ZDUM1,ZDUM2
DOUBLE PRECISION CABS1
DOUBLE PRECISION DREAL,DIMAG
COMPLEX*16 ZDUMR,ZDUMI
DREAL(ZDUMR) = ZDUMR
DIMAG(ZDUMI) = (0.0D0,-1.0D0)*ZDUMI
CSIGN(ZDUM1,ZDUM2) = CDABS(ZDUM1)*(ZDUM2/CDABS(ZDUM2))
CABS1(ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM))
PL = 1
PU = 0
IF (JOB .EQ. 0) GO TO 60
DO 20 J = 1, P
SWAPJ = JPVT(J) .GT. 0
NEGJ = JPVT(J) .LT. 0
JPVT(J) = J
IF (NEGJ) JPVT(J) = -J
IF (.NOT.SWAPJ) GO TO 10
IF (J .NE. PL) CALL ZSWAP(N,X(1,PL),1,X(1,J),1)
JPVT(J) = JPVT(PL)
JPVT(PL) = J
PL = PL + 1
10 CONTINUE
20 CONTINUE
PU = P
DO 50 JJ = 1, P
J = P - JJ + 1
IF (JPVT(J) .GE. 0) GO TO 40
JPVT(J) = -JPVT(J)
IF (J .EQ. PU) GO TO 30
CALL ZSWAP(N,X(1,PU),1,X(1,J),1)
JP = JPVT(PU)
JPVT(PU) = JPVT(J)
JPVT(J) = JP
30 CONTINUE
PU = PU - 1
40 CONTINUE
50 CONTINUE
60 CONTINUE
IF (PU .LT. PL) GO TO 80
DO 70 J = PL, PU
QRAUX(J) = DCMPLX(DZNRM2(N,X(1,J),1),0.0D0)
WORK(J) = QRAUX(J)
70 CONTINUE
80 CONTINUE
LUP = MIN0(N,P)
DO 200 L = 1, LUP
IF (L .LT. PL .OR. L .GE. PU) GO TO 120
MAXNRM = 0.0D0
MAXJ = L
DO 100 J = L, PU
IF (DREAL(QRAUX(J)) .LE. MAXNRM) GO TO 90
MAXNRM = DREAL(QRAUX(J))
MAXJ = J
90 CONTINUE
100 CONTINUE
IF (MAXJ .EQ. L) GO TO 110
CALL ZSWAP(N,X(1,L),1,X(1,MAXJ),1)
QRAUX(MAXJ) = QRAUX(L)
WORK(MAXJ) = WORK(L)
JP = JPVT(MAXJ)
JPVT(MAXJ) = JPVT(L)
JPVT(L) = JP
110 CONTINUE
120 CONTINUE
QRAUX(L) = (0.0D0,0.0D0)
IF (L .EQ. N) GO TO 190
NRMXL = DCMPLX(DZNRM2(N-L+1,X(L,L),1),0.0D0)
IF (CABS1(NRMXL) .EQ. 0.0D0) GO TO 180
IF (CABS1(X(L,L)) .NE. 0.0D0)
* NRMXL = CSIGN(NRMXL,X(L,L))
CALL ZSCAL(N-L+1,(1.0D0,0.0D0)/NRMXL,X(L,L),1)
X(L,L) = (1.0D0,0.0D0) + X(L,L)
LP1 = L + 1
IF (P .LT. LP1) GO TO 170
DO 160 J = LP1, P
T = -ZDOTC(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
CALL ZAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
IF (J .LT. PL .OR. J .GT. PU) GO TO 150
IF (CABS1(QRAUX(J)) .EQ. 0.0D0) GO TO 150
TT = 1.0D0 - (CDABS(X(L,J))/DREAL(QRAUX(J)))**2
TT = DMAX1(TT,0.0D0)
T = DCMPLX(TT,0.0D0)
TT = 1.0D0
* + 0.05D0*TT
* *(DREAL(QRAUX(J))/DREAL(WORK(J)))**2
IF (TT .EQ. 1.0D0) GO TO 130
QRAUX(J) = QRAUX(J)*CDSQRT(T)
GO TO 140
130 CONTINUE
QRAUX(J) = DCMPLX(DZNRM2(N-L,X(L+1,J),1),0.0D0)
WORK(J) = QRAUX(J)
140 CONTINUE
150 CONTINUE
160 CONTINUE
170 CONTINUE
QRAUX(L) = X(L,L)
X(L,L) = -NRMXL
180 CONTINUE
190 CONTINUE
200 CONTINUE
RETURN
END
SUBROUTINE ZQRSL(X,LDX,N,K,QRAUX,Y,QY,QTY,B,RSD,XB,JOB,INFO)
INTEGER LDX,N,K,JOB,INFO
COMPLEX*16 X(LDX,1),QRAUX(1),Y(1),QY(1),QTY(1),B(1),RSD(1),XB(1)
INTEGER I,J,JJ,JU,KP1
COMPLEX*16 ZDOTC,T,TEMP
LOGICAL CB,CQY,CQTY,CR,CXB
COMPLEX*16 ZDUM
DOUBLE PRECISION CABS1
DOUBLE PRECISION DREAL,DIMAG
COMPLEX*16 ZDUMR,ZDUMI
DREAL(ZDUMR) = ZDUMR
DIMAG(ZDUMI) = (0.0D0,-1.0D0)*ZDUMI
CABS1(ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM))
INFO = 0
CQY = JOB/10000 .NE. 0
CQTY = MOD(JOB,10000) .NE. 0
CB = MOD(JOB,1000)/100 .NE. 0
CR = MOD(JOB,100)/10 .NE. 0
CXB = MOD(JOB,10) .NE. 0
JU = MIN0(K,N-1)
IF (JU .NE. 0) GO TO 40
IF (CQY) QY(1) = Y(1)
IF (CQTY) QTY(1) = Y(1)
IF (CXB) XB(1) = Y(1)
IF (.NOT.CB) GO TO 30
IF (CABS1(X(1,1)) .NE. 0.0D0) GO TO 10
INFO = 1
GO TO 20
10 CONTINUE
B(1) = Y(1)/X(1,1)
20 CONTINUE
30 CONTINUE
IF (CR) RSD(1) = (0.0D0,0.0D0)
GO TO 250
40 CONTINUE
IF (CQY) CALL ZCOPY(N,Y,1,QY,1)
IF (CQTY) CALL ZCOPY(N,Y,1,QTY,1)
IF (.NOT.CQY) GO TO 70
DO 60 JJ = 1, JU
J = JU - JJ + 1
IF (CABS1(QRAUX(J)) .EQ. 0.0D0) GO TO 50
TEMP = X(J,J)
X(J,J) = QRAUX(J)
T = -ZDOTC(N-J+1,X(J,J),1,QY(J),1)/X(J,J)
CALL ZAXPY(N-J+1,T,X(J,J),1,QY(J),1)
X(J,J) = TEMP
50 CONTINUE
60 CONTINUE
70 CONTINUE
IF (.NOT.CQTY) GO TO 100
DO 90 J = 1, JU
IF (CABS1(QRAUX(J)) .EQ. 0.0D0) GO TO 80
TEMP = X(J,J)
X(J,J) = QRAUX(J)
T = -ZDOTC(N-J+1,X(J,J),1,QTY(J),1)/X(J,J)
CALL ZAXPY(N-J+1,T,X(J,J),1,QTY(J),1)
X(J,J) = TEMP
80 CONTINUE
90 CONTINUE
100 CONTINUE
IF (CB) CALL ZCOPY(K,QTY,1,B,1)
KP1 = K + 1
IF (CXB) CALL ZCOPY(K,QTY,1,XB,1)
IF (CR .AND. K .LT. N) CALL ZCOPY(N-K,QTY(KP1),1,RSD(KP1),1)
IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 120
DO 110 I = KP1, N
XB(I) = (0.0D0,0.0D0)
110 CONTINUE
120 CONTINUE
IF (.NOT.CR) GO TO 140
DO 130 I = 1, K
RSD(I) = (0.0D0,0.0D0)
130 CONTINUE
140 CONTINUE
IF (.NOT.CB) GO TO 190
DO 170 JJ = 1, K
J = K - JJ + 1
IF (CABS1(X(J,J)) .NE. 0.0D0) GO TO 150
INFO = J
GO TO 180
150 CONTINUE
B(J) = B(J)/X(J,J)
IF (J .EQ. 1) GO TO 160
T = -B(J)
CALL ZAXPY(J-1,T,X(1,J),1,B,1)
160 CONTINUE
170 CONTINUE
180 CONTINUE
190 CONTINUE
IF (.NOT.CR .AND. .NOT.CXB) GO TO 240
DO 230 JJ = 1, JU
J = JU - JJ + 1
IF (CABS1(QRAUX(J)) .EQ. 0.0D0) GO TO 220
TEMP = X(J,J)
X(J,J) = QRAUX(J)
IF (.NOT.CR) GO TO 200
T = -ZDOTC(N-J+1,X(J,J),1,RSD(J),1)/X(J,J)
CALL ZAXPY(N-J+1,T,X(J,J),1,RSD(J),1)
200 CONTINUE
IF (.NOT.CXB) GO TO 210
T = -ZDOTC(N-J+1,X(J,J),1,XB(J),1)/X(J,J)
CALL ZAXPY(N-J+1,T,X(J,J),1,XB(J),1)
210 CONTINUE
X(J,J) = TEMP
220 CONTINUE
230 CONTINUE
240 CONTINUE
250 CONTINUE
RETURN
END
SUBROUTINE ZAXPY(N,ZA,ZX,INCX,ZY,INCY)
COMPLEX*16 ZX(1),ZY(1),ZA
DOUBLE PRECISION DCABS1
IF(N.LE.0)RETURN
IF (DCABS1(ZA) .EQ. 0.0D0) RETURN
IF (INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
ZY(IY) = ZY(IY) + ZA*ZX(IX)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
20 DO 30 I = 1,N
ZY(I) = ZY(I) + ZA*ZX(I)
30 CONTINUE
RETURN
END
SUBROUTINE ZCOPY(N,ZX,INCX,ZY,INCY)
COMPLEX*16 ZX(1),ZY(1)
INTEGER I,INCX,INCY,IX,IY,N
IF(N.LE.0)RETURN
IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
ZY(IY) = ZX(IX)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
20 DO 30 I = 1,N
ZY(I) = ZX(I)
30 CONTINUE
RETURN
END
SUBROUTINE ZSCAL(N,ZA,ZX,INCX)
COMPLEX*16 ZA,ZX(1)
IF(N.LE.0)RETURN
IF(INCX.EQ.1)GO TO 20
IX = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
DO 10 I = 1,N
ZX(IX) = ZA*ZX(IX)
IX = IX + INCX
10 CONTINUE
RETURN
20 DO 30 I = 1,N
ZX(I) = ZA*ZX(I)
30 CONTINUE
RETURN
END
SUBROUTINE ZSWAP (N,ZX,INCX,ZY,INCY)
COMPLEX*16 ZX(1),ZY(1),ZTEMP
IF(N.LE.0)RETURN
IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
ZTEMP = ZX(IX)
ZX(IX) = ZY(IY)
ZY(IY) = ZTEMP
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
RETURN
20 DO 30 I = 1,N
ZTEMP = ZX(I)
ZX(I) = ZY(I)
ZY(I) = ZTEMP
30 CONTINUE
RETURN
END
DOUBLE PRECISION FUNCTION DZNRM2( N, ZX, INCX)
LOGICAL IMAG, SCALE
INTEGER NEXT
DOUBLE PRECISION CUTLO, CUTHI, HITEST, SUM, XMAX, ABSX, ZERO, ONE
COMPLEX*16 ZX(1)
DOUBLE PRECISION DREAL,DIMAG
COMPLEX*16 ZDUMR,ZDUMI
DREAL(ZDUMR) = ZDUMR
DIMAG(ZDUMI) = (0.0D0,-1.0D0)*ZDUMI
DATA ZERO, ONE /0.0D0, 1.0D0/
DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 /
IF(N .GT. 0) GO TO 10
DZNRM2 = ZERO
GO TO 300
10 ASSIGN 30 TO NEXT
SUM = ZERO
NN = N * INCX
DO 210 I=1,NN,INCX
ABSX = DABS(DREAL(ZX(I)))
IMAG = .FALSE.
GO TO NEXT,(30, 50, 70, 90, 110)
30 IF( ABSX .GT. CUTLO) GO TO 85
ASSIGN 50 TO NEXT
SCALE = .FALSE.
50 IF( ABSX .EQ. ZERO) GO TO 200
IF( ABSX .GT. CUTLO) GO TO 85
ASSIGN 70 TO NEXT
GO TO 105
100 ASSIGN 110 TO NEXT
SUM = (SUM / ABSX) / ABSX
105 SCALE = .TRUE.
XMAX = ABSX
GO TO 115
70 IF( ABSX .GT. CUTLO ) GO TO 75
110 IF( ABSX .LE. XMAX ) GO TO 115
SUM = ONE + SUM * (XMAX / ABSX)**2
XMAX = ABSX
GO TO 200
115 SUM = SUM + (ABSX/XMAX)**2
GO TO 200
75 SUM = (SUM * XMAX) * XMAX
85 ASSIGN 90 TO NEXT
SCALE = .FALSE.
HITEST = CUTHI/FLOAT( N )
90 IF(ABSX .GE. HITEST) GO TO 100
SUM = SUM + ABSX**2
200 CONTINUE
IF(IMAG) GO TO 210
ABSX = DABS(DIMAG(ZX(I)))
IMAG = .TRUE.
GO TO NEXT,( 50, 70, 90, 110 )
210 CONTINUE
DZNRM2 = DSQRT(SUM)
IF(SCALE) DZNRM2 = DZNRM2 * XMAX
300 CONTINUE
RETURN
END
COMPLEX*16 FUNCTION ZDOTC(N,ZX,INCX,ZY,INCY)
COMPLEX*16 ZX(1),ZY(1),ZTEMP
ZTEMP = (0.0D0,0.0D0)
ZDOTC = (0.0D0,0.0D0)
IF(N.LE.0)RETURN
IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
IX = 1
IY = 1
IF(INCX.LT.0)IX = (-N+1)*INCX + 1
IF(INCY.LT.0)IY = (-N+1)*INCY + 1
DO 10 I = 1,N
ZTEMP = ZTEMP + DCONJG(ZX(IX))*ZY(IY)
IX = IX + INCX
IY = IY + INCY
10 CONTINUE
ZDOTC = ZTEMP
RETURN
20 DO 30 I = 1,N
ZTEMP = ZTEMP + DCONJG(ZX(I))*ZY(I)
30 CONTINUE
ZDOTC = ZTEMP
RETURN
END
DOUBLE PRECISION FUNCTION DCABS1(Z)
COMPLEX*16 Z,ZZ
DOUBLE PRECISION T(2)
EQUIVALENCE (ZZ,T(1))
ZZ = Z
DCABS1 = DABS(T(1)) + DABS(T(2))
RETURN
END
Done. The code is compiled via
f2py. We even get a nice __doc__ generated.
In [5]:
print(mtm_mann.__doc__)
The next cell defines a python wrapper (to call the code we just compiled as a module) and a function to plot the results.
In [6]:
def compute_mtm(arr, dt=1.0, npi=2, nwin=3, f1=0.0, f2=0.0, inorm=0, ispec=1,
iresh=1, ithresh=3, inoise=0, ilog=1, ismooth=1, isignal=0,
irecon=0, nsignals=0, iplotresh=1, iplotftest=1, iplotsmoo=1,
iplotraw=1, iplotconf=1, icon=0):
if nwin > (2 * npi - 1):
nwin = 2 * npi - 1
if f2 == 0:
f2 = 0.5 / dt
if isignal == 2:
inoise = 2
iplotsmoo = 0
iplotconf = 0
if isignal == 1:
iplotftest = 0
iplotresh = 0
spec_raw, spec_resh, spec_smoo, spec_conf, recon_sig = mtm_mann(
arr, dt, npi, nwin, f1, f2, inorm, ispec, iresh, ithresh, inoise, ilog,
ismooth, isignal, irecon, nsignals, iplotresh, iplotftest, iplotsmoo,
iplotraw, iplotconf, icon)
spec_raw = resize_spec(spec_raw)
spec_resh = resize_spec(spec_resh)
spec_smoo = resize_spec(spec_smoo)
spec_conf = resize_spec(spec_conf)
recon_sig = resize_spec(recon_sig)
return spec_raw, spec_resh, spec_smoo, spec_conf, recon_sig
def plot_spec(spec_raw):
fig, ax = plt.subplots()
m = np.array([1.0, 2.0, 5.0, 10.0, 20.0, 100.0])
ohm = 1.0 / m
y = 110
ax.loglog(spec_raw[:, 0], spec_raw[:, 1], linewidth=0.5)
for i in range(0, 6):
ax.text(ohm[i], y, str(int(m[i])))
ax.set_xlabel(r'Frequency [years $^{-1}$]')
ax.set_ylabel('Spectral Power')
ax.set_title('MTM Spectrum\n',)
ax.text(2, y, 'Period [years]')
ax.text(.01, 0.01, r'Nino3 from Speedy')
ax.grid(True)
return fig, ax
def resize_spec(arr):
arr = np.asanyarray(arr)
if arr.ndim != 2:
raise ValueError("Array must be 2D.")
size = 1
nrow, ncol = arr.shape
for i in range(1, nrow):
if arr[i, 0] != 0:
size = size + 1
if size != 1:
arr = arr[0:size, 0:ncol]
return arr
Let's test the module.
In [7]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
X = [1.0135518, -0.7113242, -0.3906069, 1.565203, 0.0439317, -1.1656093,
1.0701692, 1.0825379, -1.2239744, -0.0321446, 1.1815997, -1.4969448,
-0.7455299, 1.0973884, -0.2188716, -1.0719573, 0.9922009, 0.4374216,
-1.6880219, 0.2609807]
X = np.atleast_2d(X)
spec_raw, spec_resh, spec_smoo, spec_conf, recon_sig = compute_mtm(X)
fig, ax = plot_spec(spec_raw)
When compiling code on the fly like this we can easily debug, develop, test and plot the results. All in a literate notebook! Now, from Futurama, a well deserved beer:
In [8]:
from IPython.display import Image
Image("http://s18.postimg.org/xoghfy291/Olde_Fortran_Malt_Liquor.jpg")
Out[8]:
(Simpsons and Futurama writers must be Fortran programers!)
In [9]:
HTML(html)
Out[9]:
This post was written as an IPython notebook. It is available for download or as a static html.
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.