c     Programme beta_eca.f
c     *********************
c
c     Ce programme compare deux sources photometriques dans
c     " beta.mes ".
c     Ce programme supprime certaines etoiles de la compa-
c     raison.
c     Ce programme fournit la moyenne et l'ecart-type des
c     differences de beta entre les deux references etudiees.
c
c
c     Sortie dans les fichiers : " beta.eca " et "graph.eca".
c
c
c
      character stat*3
      integer etoil, eto, ngde(16)
      logical l
      dimension num(150),beta1(150),beta2(150),delta(150)
c
      data ne/0/,etoil/1/,n/1/
      data ngde/0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/
c
c
      open(10,file='.T/beta.mes',status='old')
      inquire(file='.T/beta.eca',exist=l)
      if(.not.l) then
        stat='new'
        else
        stat='old'
      endif
      open(15,file='.T/beta.eca',status=stat)
c
c
c     Sortie graphique : ".T/graph.eca"
c
      inquire(file='.T/graph.eca',exist=l)
      if(.not.l)then
        stat='new'
        else
        stat='old'
      endif
      open(20,file='.T/graph.eca',status=stat)
c
c     lecture des donnees
c
      print 300
  300 format(/,'Entrez la premiere reference: ',$)
      read *,iref1
      print 301
  301 format(/,'Entrez la deuxieme reference: ',$)
      read *,iref2
      print 302
  302 format(/,'Donnez le nombre d''etoiles a supprimer: ',$)
      read *,nreb
      if(nreb.eq.0) goto 2
      print 303
  303 format(/,'Donnez le numero des etoiles: ')
      read *,(ngde(j),j=1,nreb)
c
c
    2 eto=etoil
    3 read(10,100,end=11) etoil,iref,beta
  100 format(1x,i4,1x,i4,1x,f6.3)
c   
      do 35 j=1,nreb
         if(etoil.eq.ngde(j)) goto 3
   35 continue
c
    4 if(eto-etoil) 10,6,10
c
    6 if(ne.eq.1) goto 8
      if(ne.eq.0) goto 7
c
    7 if(iref.ne.iref1) goto 3
      ne=1
      bet1=beta
      goto 3
c
    8 if(iref.ne.iref2) goto3
      bet2=beta
      beta1(n)=bet1
      beta2(n)=bet2
      num(n)=etoil
      delta(n)=bet1-bet2
      n=n+1
      ne=0
      goto 3
c
   10 ne=0
      eto=etoil
      goto 4
c
c     Sortie de la premiere partie
   11 maximum=abs(n-1)
c
c
      write(15,109)
      write(15,110)
  109 format(/,'    No     ref1    ref2   delta ')
  110 format('  -------------------------------')
      do 89 j=1,maximum
      write(15,101)num(j),beta1(j),beta2(j),delta(j)
      write(20,101)num(j),beta1(j),beta2(j),delta(j)
   89 continue
  101 format(2x,i4,3x,f6.3,2x,f6.3,2x,f6.3)
c
c
c     Premier passage
c
      deltam=0
      do 36 j=1,maximum
         deltam=deltam+delta(j)
   36 continue
c
      deltam=deltam/maximum
      sigm=0
      do 37 j=1,maximum
         sigm=(delta(j)-deltam)**2+sigm
   37 continue
      sigm=SQRT(sigm/(maximum-1))
c
c
      print 103
  103 format(/,3x,'Premier passage : ')
      print 111
  111 format(3x,'-----------------')
      print 102,deltam,sigm
  102 format(/,3x,'delta moyen = ',f7.4,/,
     &         3x,'sigma =       ',f7.4)
      print 108,maximum
  108 format(/,3x,'Nombre d''etoiles considerees : ',i4)
      if(sigm.LT.0.02) goto 121
c
c
c     Deuxieme passage
c
      n2=0
      deltam2=0
      do 38 j=1,maximum
         if(abs(delta(j)).GT.(2*sigm)) goto 39
         deltam2=deltam2+delta(j)
         n2=n2+1
   39 continue
   38 continue
      maximum2=n2
      deltam2=deltam2/maximum2
c     sigm2=0
      sigm2=0
      do 40 j=1,maximum
         if(abs(delta(j)).GT.(2*sigm)) goto 41
         sigm2=sigm2+(delta(j)-deltam2)**2
   41 continue
   40 continue
      sigm2=SQRT(sigm2/(maximum2-1))
c
      print 104
  104 format(/,3x,'Deuxieme passage : ')
      print 112
  112 format(3x,'------------------')
      print 102,deltam2,sigm2
      print 105,maximum2
  105 format(3x,'Nombre d''etoiles considerees :  ',i3)
c
c
c     Troisieme passage
c
      if(sigm.GT.0.03.and.(sigm/sigm2).GT.1.5) then
c
      n3=0
      deltam3=0
      sigm3=0
c
      do 42 j=1,maximum
         if(abs(delta(j)).GT.(2*sigm2)) goto 43
         deltam3=deltam3+delta(j)
         n3=n3+1
   43 continue
   42 continue
      maximum3=n3
      deltam3=deltam3/maximum3
c
      do 44 j=1,maximum
         if(abs(delta(j)).GT.(2*sigm2)) goto 45
         sigm3=sigm3+(delta(j)-deltam3)**2
   45 continue
   44 continue
      sigm3=SQRT(sigm3/(maximum3-1))
c
      print 106
  106 format(/,3x,'Troisieme passage : ')
      print 113
  113 format(3x,'------------------')
      print 102,deltam3,sigm3
      print 107,maximum3
  107 format(3x,'Nombre d''etoiles considerees : ',i3)
c
      endif
c
  121 stop
c
      end