← Back to branch summary

~ubuntu-branches/debian/sid/octave-tisean/sid

« back to all changes in this revision

Viewing changes to src/source_f/store_spec.f

  • Committer: Package Import Robot
  • Author(s): Rafael Laboissiere
  • Date: 2017-08-14 12:53:47 UTC
  • Revision ID: package-import@ubuntu.com-20170814125347-ju5owr4dggr53a2n
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/source_f/store_spec.f
 
1
c===========================================================================
 
2
c
 
3
c   This file is part of TISEAN
 
4
c 
 
5
c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
 
6
c 
 
7
c   TISEAN is free software; you can redistribute it and/or modify
 
8
c   it under the terms of the GNU General Public License as published by
 
9
c   the Free Software Foundation; either version 2 of the License, or
 
10
c   (at your option) any later version.
 
11
c
 
12
c   TISEAN is distributed in the hope that it will be useful,
 
13
c   but WITHOUT ANY WARRANTY; without even the implied warranty of
 
14
c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
15
c   GNU General Public License for more details.
 
16
c
 
17
c   You should have received a copy of the GNU General Public License
 
18
c   along with TISEAN; if not, write to the Free Software
 
19
c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 
20
c
 
21
c===========================================================================
 
22
c   store data periodogram of x 
 
23
c   if iback.ne.0 transform back to get autocorrelation instead
 
24
C   author Thomas Schreiber (1998)
 
25
c===========================================================================
 
26
      subroutine store_spec(nmax,x,iback)
 
27
      parameter(nx=1000000)
 
28
      dimension x(nmax), w1(nx), w2(nx), iw(15)
 
29
      save w2, iw
 
30
 
 
31
      if(nmax.gt.nx) then
 
32
         call xstopx ("store_spec: make nx larger.")
 
33
      endif
 
34
      call rffti1(nmax,w2,iw)  
 
35
      call rfftf1(nmax,x,w1,w2,iw)
 
36
      do 10 n=1,nmax
 
37
 10      x(n)=x(n)/real(nmax)
 
38
      x(1)=x(1)**2
 
39
      do 20 n=2,(nmax+1)/2
 
40
         amp=x(2*n-2)**2+x(2*n-1)**2
 
41
         pha=atan2(x(2*n-1),x(2*n-2))
 
42
         x(2*n-2)=amp
 
43
 20      x(2*n-1)=pha
 
44
      if(mod(nmax,2).eq.0) x(nmax)=x(nmax)**2
 
45
      if(iback.eq.0) return
 
46
      do 30 n=1,nmax
 
47
 30      x(n)=x(n)*nmax
 
48
      do 40 n=2,(nmax+1)/2
 
49
 40      x(2*n-1)=0
 
50
      call rfftb1(nmax,x,w1,w2,iw)
 
51
      end