C PROGRAMA spectra1.for *************************************
C Tabare Gallardo, abril/98, gallardo@fisica.edu.uy
C Facultad de Ciencias - Obs. Astr. Los Molinos
C Calcula espectro de frecuencias en un cierto intervalo
C FUNCION MODELO  Y = C1 + C2*COS(2PIwt) + C3*SIN(2PIwt)
C Pesos en forma Hanning para eliminar lobulos secundarios

      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION T(10000),X(10000),P(10000),xp(10000)
      COMMON T,XP,P,N,TWOPI,F1X
      TWOPI=DATAN(1.D0)*8.D0
          
      write(*,*) 'ESPECTRO DE LA SERIE TEMPORAL datos.dat ', 
     *'(equiespaciada o no)'   
      write(*,*) 'AJUSTANDO AL MODELO Y = C1 + C2*COS(2PI*f*t) + ',
     *'C3*SIN(2PI*f*t)'
      write(*,*) 'Refs: Astronomical Journal (1981) vol.86(4), pp 619.'
      write(*,*) '      Astronomical Journal (1996) vol.111(1), pp 541.'
      write(*,*) '      Astronomical Journal (1997) vol.113(2), pp 863.'
      write(*,*)      
      write(*,*)'El espectro estara en spectra1.dat y el orden sera:'
      write(*,*) 'frecuencia, Coef.Corr.Espec., C1, C2, C3'
      write(*,*)       
      write(*,*) 'Nota: unidad de frecuencia = ciclos por ', 
     *'unidad de tiempo'
      write(*,*)
 
      OPEN(1,FILE='datos.dat',STATUS='OLD')
      do 50 i=1,10000
        read(1,*,end=52) t(i),x(i)
  50  continue  
  52  close(1)
      n=i-1           
      write(*,*) 'son ',n,'    datos en la serie temporal'   
      
      zn=t(n)-t(1)
      ancho=2.d0/zn 
      write(*,*) 'ancho de cada linea espectral =',ancho  
      write(*,*)
      write(*,*) 'Definicion de los limites del espectro:'
      write(*,*) 'frecuencia inicial?'
      read(*,*) frecini      
      write(*,*) 'frecuencia final?'
      read(*,*) frecfin  
      write(*,*) 'Numero de puntos que tendra el espectro?' 
      read(*,*) nwf
      dnwf1=dfloat(nwf-1)
      FFFF=(frecfin-frecini)/dnwf1          
      tau=t(1)+zn/2.d0
      
C DEFINO PESOS CENTRADO EN TAU
      sumapeso=0.d0
      DO 10 I=1,N
        ttau=t(i)-tau
c hanning window
        p(i)=1.d0+dcos(twopi*ttau/zn)
        sumapeso=sumapeso+p(i)
  10  CONTINUE

c normalizo pesos y defino XP, F1X y varianza
      F1X=0.d0         
      xcuad=0.d0
      do 13 i=1,n
        p(i)=p(i)/sumapeso
        xp(i)=x(i)*p(i)   
        f1x=f1x+xp(i)    
        xcuad=xcuad+x(i)*xp(i)
  13  continue
      varianza=xcuad-f1x*f1x
                                                
      write(*,*) 'estoy calculando el espectro....'                   
                 
      ancho100=ancho/1.d2                          
      open(2,file='spectra1.dat',status='unknown')         
      do 41 nw=1,nwf
        W=frecini+FFFF*dfloat(nw-1)   
c evito calculo de frec proxima a cero
          if(dabs(w).lt.ancho100) then
             goto 41
          end if
        call fosyl(W,C1,C2,C3,YY)
C AMPLITUD TEORICA de la sinusoide = DSQRT(C2*C2+C3*C3)
C CENTRO O NIVEL CERO = C1
C coeficiente de correlacion espectral (Sylvio-Foster) = CCE    
        CCE=(YY-F1X*F1X)/varianza
        write(2,11) W,CCE,C1,C2,C3
  41  continue
  11  FORMAT(e14.7,e16.9,3e14.7)
      close(2)           
      END

C SUBRUTINA FOSYL (FOSTER-SYLVIO)********************************
C FUNCION MODELO Y = C1 + C2*COS(2PIwt) + C3*SIN(2PIwt)
      SUBROUTINE FOSYL(W,C1,C2,C3,YY)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION XP(10000),T(10000),P(10000)
      DIMENSION CA(10000),SA(10000)
      COMMON T,XP,P,N,TWOPI,F1X
      F11=1.D0
      F12=0.D0
      F13=0.D0
      F22=0.D0
      F23=0.D0
      F33=0.D0
      F2X=0.D0
      F3X=0.D0
      DO 20 I=1,N
      ALFA=TWOPI*T(I)*W
      AL=DMOD(ALFA,TWOPI)
      CA(I)=DCOS(AL)
      SA(I)=DSIN(AL)
      CAP=CA(I)*P(I)
      SAP=SA(I)*P(I)
      F12=F12+CAP
      F13=F13+SAP
      F22=F22+CA(I)*CAP
      F23=F23+CA(I)*SAP
      F33=F33+SA(I)*SAP
C      F1X=F1X+XP(I)      
C F1X es igual al valor medio ponderado de X y de Y
      F2X=F2X+XP(I)*CA(I)
      F3X=F3X+XP(I)*SA(I)
   20 CONTINUE
      DET=F11*(F22*F33-F23*F23)
      DET=DET-F12*(F12*F33-F23*F13)
      DET=DET+F13*(F12*F23-F22*F13)
      C1=-F12*(F2X*F33-F23*F3X)+F13*(F2X*F23-F22*F3X)
      C1=(C1+F1X*(F22*F33-F23*F23))/DET
      C2=F11*(F2X*F33-F23*F3X)+F13*(F12*F3X-F2X*F13)
      C2=(C2-F1X*(F12*F33-F23*F13))/DET
      C3=F11*(F22*F3X-F2X*F23)-F12*(F12*F3X-F2X*F13)
      C3=(C3+F1X*(F12*F23-F22*F13))/DET
C CALCULO YY
      YY=0.D0
      DO 21 I=1,N
      Y=C1+C2*CA(I)+C3*SA(I)
      YY=YY+Y*Y*P(I)
   21 CONTINUE
      RETURN
      END      
C *********************************************************