C @(#)tsaaov.for 17.1.1.1 (ES0-DMD) 01/25/02 17:20:36 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C @(#)tsaaov.for 5.1 (ESO-IPG) 4/5/93 15:58:45 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT (c) 1992 European Southern Observatory and Warsaw Observatory C.IDENT tsaaov.for C.AUTHOR Alex Schwarzenberg-Czerny, ESO and Warsawe Observatory C.KEYWORD MIDAS, time series, AOV/TSA C.LANGUAGE FORTRAN 77 C.PURPOSE Compute analysis of variance periodogramme C.VERSION 0.0 June 1992 C.RETURNS None C.ENVIRON TIME/ANAL context C----------------------------------------------------------------------------- C C INCLUDE 'MID_REL_INCL:TSA_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DEF.INC' C INTEGER NAXIS PARAMETER (NAXIS=2) REAL*8 START(NAXIS) ! START FREQUENCY REAL*8 STEP(NAXIS) ! STEP IN FREQUENCY INTEGER LENGTH(NAXIS) ! NUMBER OF FREQUENCIES INTEGER ORDER ! NUMBER OF BINS INTEGER COVER ! NUMBER OF BIN COVERAGES CHARACTER*60 INAME ! NAME OF INPUT TABLE CHARACTER*60 ONAME ! NAME OF OUTPUT IMAGE C INTEGER IACTS,KUN,KNUL INTEGER TID,ITIME,IVAL INTEGER NCOL,ICOL,NROW,IROW,ISOR INTEGER LFIELD,TTYP,VTYP INTEGER PTIME,PVAL,PPER,IDPER CHARACTER*10 FORM CHARACTER*80 TEXT C INCLUDE 'MID_REL_INCL:TSA_DAT.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C Get parameters C CALL STSPRO ('tsaaov') CALL STKRDC ('IN_A' ,1,1,60,IACTS,INAME ,KUN,KNUL,ISTAT) CALL STKRDC ('OUT_A', 1,1,60,IACTS,ONAME ,KUN,KNUL,ISTAT) CALL STKRDD ('STARTTSA', $ 1, 1,IACTS,START(1) ,KUN,KNUL,ISTAT) CALL STKRDD ('STEPTSA', $ 1, 1,IACTS,STEP(1) ,KUN,KNUL,ISTAT) CALL STKRDI ('NSTEPS', 1, 1,IACTS,LENGTH(1),KUN,KNUL,ISTAT) CALL STKRDI ('ORDERTSA', $ 1, 1,IACTS,ORDER, KUN,KNUL,ISTAT) CALL STKRDI ('COVER', 1, 1,IACTS,COVER, KUN,KNUL,ISTAT) START(NAXIS)= ZERO STEP(NAXIS)= ONE LENGTH(NAXIS)= 2 C C Map input data C CALL TBTOPN (INAME,F_I_MODE,TID,ISTAT) CALL TBIGET (TID,NCOL,NROW,ISOR,ICOL,IROW,ISTAT) CALL TBLSER (TID,'TIME',ITIME,ISTAT) IF (ITIME.LT.0) THEN CALL STETER(3,'Column :TIME not found') ENDIF CALL TBLSER (TID,'VALUE',IVAL,ISTAT) IF (IVAL.LT.0) THEN CALL STETER(3,'Column :VALUE not found') ENDIF CALL TBFGET (TID,ITIME,FORM,LFIELD,TTYP,ISTAT) CALL TBFGET (TID,IVAL, FORM,LFIELD,VTYP,ISTAT) CALL TBDGET (TID,ISTORE,ISTAT) IF (ISTORE.NE.F_TRANS) THEN TEXT='Input table '//INAME//' stored not transposed' CALL STETER(2,TEXT) ENDIF IF (TTYP.NE.D_R8_FORMAT.OR.VTYP.NE.D_R8_FORMAT) THEN CALL STETER(1,'Column(s) must be of DOUBLE PRECISION type') ENDIF CALL TBCMAP (TID,ITIME,PTIME,ISTAT) CALL TBCMAP (TID,IVAL, PVAL, ISTAT) C C Map output data C CALL STIPUT (ONAME,D_R8_FORMAT,F_IO_MODE, $ F_IMA_TYPE,NAXIS,LENGTH, $ START,STEP,ONAME,'AXIS: 1/TIME DATA: UNITLESS', $ PPER,IDPER,ISTAT) C C Compute AOV periodogramme C CALL TIMAOV (MADRID(PTIME),MADRID(PVAL), $ MADRID(PPER),NROW,LENGTH(1), $ START(1),STEP(1),ORDER,COVER) C C Wind-up C CALL DSCUPT(IDPER,IDPER,' ',ISTAT) CALL STSEPI C END C C C SUBROUTINE TIMAOV(T,X,Y,N,NFREQ,FREQIN,FSTEP,NORDER,NCOVER) C C Analysis of Variance Periodogramme C Reference: M.N.R.A.S. 241, 153 C Alex Schwarzenberg-Czerny Warsaw Observatory 1988 C C INCLUDE 'MID_REL_INCL:TSA_DEF.INC' C INTEGER NORDER ! ORDER (NUMBER OF BINS) INTEGER NCOVER ! NUMBER OF BIN COVERS INTEGER N ! NUMBER OF OBSERVATIONS INTEGER NFREQ ! NUMBER OF FREQUENCIES DOUBLE PRECISION T(N),X(N) ! TIMES AND VALUES OF OBSERVATIONS DOUBLE PRECISION Y(NFREQ,2) ! PERIODOGRAMME REAL*8 FREQIN,FSTEP ! FREQUENCY START AND INCREMENT C INTEGER MAXORD PARAMETER(MAXORD=100) C INTEGER ICOVER,IF,I,NO,IBIN,NR INTEGER NINCMPL,NINSUFF,NALL INTEGER NBIN(MAXORD) DOUBLE PRECISION BMEAN(MAXORD) DOUBLE PRECISION T0,XM,VA,SAVE,COVER,VAR DOUBLE PRECISION FREQ,XMEAN,S1,S2 CHARACTER*80 TEXT C INCLUDE 'MID_REL_INCL:TSA_CONST.INC' C IF (NORDER.GT.MAXORD) THEN WRITE(TEXT,'(A,I3)') 1 'WARNING: Used maximum alllowed order of ',MAXORD CALL STTPUT(TEXT,ISTAT) NORDER=MAXORD ENDIF C C Set variables (incl. mean and variance) C NINCMPL=0 NINSUFF=0 T0=T(1) XM=ZERO VA=ZERO DO 10 I=1,N XM=XM+X(I) 10 CONTINUE XM=XM/N DO 11 I=1,N SAVE=X(I)-XM VA=VA+SAVE*SAVE 11 CONTINUE DO 12 IF=1,NFREQ DO 14 I=1,2 Y(IF,I)=ZERO 14 CONTINUE 12 CONTINUE C C Loop over coverages and frequencies C DO 1 ICOVER=1,NCOVER COVER=(ICOVER-1.)/NCOVER NO=NORDER DO 7 IF=1,NFREQ FREQ=((IF-1)*FSTEP+FREQIN)*NO DO 3 IBIN=1,NO BMEAN(IBIN)=ZERO NBIN(IBIN)=0 3 CONTINUE C C Assign data to bins and add C DO 2 I=1,N IBIN=1+MOD(INT((T(I)-T0)*FREQ+COVER),NO) BMEAN(IBIN)=X(I)+BMEAN(IBIN) NBIN(IBIN)=1+NBIN(IBIN) 2 CONTINUE C C Correct for incomplete bins C NR=NO NALL=N XMEAN=0. VAR=VA DO 4 IBIN=1,NO BMEAN(IBIN)=BMEAN(IBIN)-NBIN(IBIN)*XM IF (NBIN(IBIN).LE.1) THEN NINCMPL=NINCMPL+1 NR=NR-1 XMEAN=XMEAN-BMEAN(IBIN) VAR=VAR-BMEAN(IBIN)*BMEAN(IBIN) NALL=NALL-NBIN(IBIN) ENDIF 4 CONTINUE C C Calculate A.O.V. statistics for a given frequency C IF (NR.GT.1) THEN XMEAN=XMEAN/NALL S1=ZERO DO 5 IBIN=1,NO IF (NBIN(IBIN).GT.1) THEN SAVE=BMEAN(IBIN)/NBIN(IBIN)-XMEAN S1=SAVE*SAVE*NBIN(IBIN)+S1 ENDIF 5 CONTINUE IF (S1.LT.VAR) THEN S2=VAR-S1 Y(IF,1)=Y(IF,1)+S1/(NR-1)/S2*(NALL-NR) ELSE NINSUFF=NINSUFF+1 Y(IF,1)=Y(IF,1)+ONE Y(IF,2)=MAX(Y(IF,2),ONE) IF (NINSUFF.LE.1) THEN CALL STTPUT( $ 'fatal calculation/rounding error',ISTAT) ENDIF ENDIF ELSE NINSUFF=NINSUFF+1 Y(IF,1)=Y(IF,1)+ONE Y(IF,2)=MAX(Y(IF,2),TWO) ENDIF 7 CONTINUE 1 CONTINUE DO 20 IF=1,NFREQ Y(IF,1)=Y(IF,1)/NCOVER 20 CONTINUE IF (NINCMPL.GT.0.OR.NINSUFF.GT.0) THEN CALL STTPUT( $ 'Statistics of bad events (see quality row for details):', $ ISTAT) WRITE(TEXT,'(2(i8,a))') $ NINCMPL,' underpopulated bins' CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,'(2(i8,a))') $ NINSUFF,' of',NFREQ*NCOVER, $ ' bad frequencies' CALL STTPUT(TEXT,ISTAT) ENDIF Y(NFREQ,2)=TWO END C