c     program bda/progf/dered.f
c
c======================================================================= 
c     computes for each star: the color excesses        e(b-v)&e(u-b)  =
c                             the intrinsic colors      (b-v)0&(u-b)0  =
c                             the apparent visual corrected mag    v0  =
c                             the intrinsic visual mag            mv0  =
c                             the distance modelus            (v0-mv0) =
c            for the cluster: the average color excesses e(b-v)&e(u-b) = 
c                             the average distance modelus   (v0-mv0)  =
c            using the relations: (u-b)0=3.69(b-v)0+0.03               =
c                                 e(u-b)/e(b-v)=0.69+0.05e(b-v)        =
c                                 av=3.2e(b-v)                         =
c                                                                      =
c     programme apporte par Aiad et developpe avec lui en avril 88     =
c     version simplifiee pour calculer les couleurs derougies          =
c     pour diagrammes UBV                              24.5.89         =
c=======================================================================
c
      dimension v(200),bmv(200),umb(200),ebmv(200),eumb(200)
      dimension v0(200),bmv0(200),umb0(200),xmv0(200),dm(200)
      dimension xmu(11),xmb(20),ub(11),bv(20),ubsg(14),bvsg(14)
      dimension difdm(200),eumb2(200),ebmv2(200),ebdif(200)
      dimension bmv02(200),umb02(200),sebmv(200),seumb(200),sdm(200)
      dimension adm(200),no(200)
      integer sn(200),ref(200),memb(200)
      integer finp
      character stat*3,alp*5,bet*5,afr*5
      logical ll
c
      data xmu/-2.12,-1.6,-1.04,-0.52,-0.03,0.34,0.75,1.02,1.26
     1,1.62,1.84/
      data ub/-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0
     1,0.05/
      data xmb/0.60,1.15,1.55,1.8,2.0,2.25,2.7,2.95,3.2,3.55,4.20
     1,4.80,5.40,5.90,6.30,6.80,7.20,7.65,8.11,8.60/
      data bv/-0.10,-0.05,0.0,0.05,0.10,0.15,0.25,0.30,0.35,0.40,0.50
     1,0.60,0.70,0.80,0.90,1.00,1.10,1.20,1.30,1.40/
      data ubsg/-1.17,-1.16,-1.14,-1.13,-1.13,-1.06,-1.00,-0.93,-0.83
     1,-0.72,-0.69,-0.63,-0.55,-0.49/
      data bvsg/-0.31,-0.31,-0.31,-0.29,-0.27,-0.23,-0.19,-0.17,-0.13
     1,-0.10,-0.08,-0.05,-0.03,-0.02/
c
      data finp/1/,fr/3.2/
c
c     parametres alpha et beta des "droites" de rougissement
c
      narg=iargc()
      if(narg.eq.0) then
        alpha=0.70
        beta=0.04
      else
        call getarg(1,alp)
        read(alp,500) alpha
        call getarg(2,bet)
        read(bet,500) beta
500     format(f5.2)
        if(narg.eq.3) then
          call getarg(3,afr)
          read(afr,500) fr
          endif
      endif
c
c     reading the source file
c
      open(10,file='ubv.cmd',status='old')
      inquire(file='tmp.drd',exist=ll)
      if(.not.ll) then
        stat='new'
      else
        stat='old'
      endif
      open(18,file='tmp.drd',status=stat)
      ipas=0
      i2=0
      im=0
      i=0
19    i=i+1
       read(10,20,end=31) no(i),ref(i),v(i),bmv(i),umb(i)
20     format(1x,i4,1x,i4,1x,f7.3,2(1x,f6.3))
       if(bmv(i).gt.1.8) then
       i2=i2+1
       goto 19
       endif
       if(umb(i).ge.0.5) goto 19
       memb(i)=1
       im=im+1
      goto 19
31    n=(i-1)-i2
      if(n.eq.0) then
        write(6,710)
710     format(/,' sorry, no stars left',/)
        stop
      endif
      if(im.eq.0) then
        write(6,700)
700     format(/,' sorry, no blue stars',/)
        stop
      endif
c
c     computation of "parameters" for all stars
c
      do 60 i=1,n
       ebmve=0.0
       ebmv(i)=ebmve
       eumb2e=0.0
       eumb2(i)=eumb2e
       k=1
       k2=1
40     eumb(i)=alpha*ebmv(i)+beta*ebmv(i)**2
       umb0(i)=umb(i)-eumb(i)
       bmv0(i)=(umb0(i)-0.03)/3.69
       if(k.gt.10) go to 50
       k=k+1
       ebmv(i)=bmv(i)-bmv0(i)
       go to 40
c
50     ebmv2(i)=1.42857*eumb2(i)
       bmv02(i)=bmv(i)-ebmv2(i)
       xx=bmv02(i)
       umb02(i)=0.0016+1.7721*xx-11.3363*xx*xx+37.1916*xx**3
     1-88.1025*xx**4+98.2920*xx**5+25.9705*xx**6-75.8215*xx**7
     2-61.7172*xx**8+78.5286*xx**9
       k2=k2+1
       if(k2.gt.10) goto 52
       eumb2(i)=umb(i)-umb02(i)
       goto 50
c
52     v0(i)=v(i)-fr*ebmv(i)
       bmv0(i)=bmv(i)-ebmv(i)
       umb0(i)=umb(i)-eumb(i)
       write(18,21) no(i),ref(i),v0(i),bmv0(i),umb0(i),ebmv(i)
21     format(1x,i4,1x,i4,1x,f7.3,3(1x,f6.3))
60     continue
c
c     fin provisoire
c
      if(finp.eq.1) stop
c
c     testing for second pass
c
61    do 620 i=1,n
        if(ipas.eq.0.or.umb0(i).ge.0.0.or.bmv0(i).ge.0.0) go to 620
        memb(i)=1
620   continue
c
c     computation of the absolute magnitude from the zams
c
      do 63 i=1,n
51     if(bmv0(i).lt.-0.1) go to 56
       do 54 l=1,20
        if(bmv0(i).gt.bv(l)) go to 54
        if(l.eq.1) then
        xmv0(i)=((xmb(l)-xmb(l+1))/(bv(l)-bv(l+1)))*
     1     (bmv0(i)-bv(l))+xmb(l)
        else
        xmv0(i)=((xmb(l)-xmb(l-1))/(bv(l)-bv(l-1)))*
     1     (bmv0(i)-bv(l-1))+xmb(l-1)
        endif
        go to 62
54    continue
c
56    do 58 l1=1,10
       if(umb0(i).gt.ub(l1).and.l1.ne.1) go to 58
       if(l1.eq.1) then
57     xmv0(i)=((xmu(l1)-xmu(l1+1))/(ub(l1)-ub(l1+1)))*
     1    (umb0(i)-ub(l1))+xmu(l1)
       else
       xmv0(i)=((xmu(l1)-xmu(l1-1))/(ub(l1)-ub(l1-1)))*
     1   (umb0(i)-ub(l1-1))+xmu(l1-1)
       endif
       go to 62
58    continue
62    dm(i)=v0(i)-xmv0(i)
63    continue
c
c     computation of mean redenning values
c
      i1=0
      iter=0
      do 90 i=1,n
       if(memb(i).ne.1) go to 90
       i1=i1+1
       sebmv(i1)=ebmv(i)
       seumb(i1)=eumb(i)
       sdm(i1)=dm(i)
90    continue
      it=i1
91    call moment(sebmv,it,ebmvm,adev1,sdev1,var1,skew1,curt1)
      call moment(seumb,it,eumbm,adev2,sdev2,var2,skew2,curt2)
      call moment(sdm,it,ddm,adev3,sdev3,var3,skew3,curt3)
      print *,' iteration: ',iter,' e(b-v) = ',ebmvm,'  it = ',it
      is=0
      do 95 i=1,n
      if(memb(i).ne.1) go to 95
      dif1=ebmv(i)-ebmvm
      dif1=abs(dif1)
      twosig=2.0*sdev1
      if(dif1.lt.twosig) then
       is=is+1
       sebmv(is)=ebmv(i)
       seumb(is)=eumb(i)
       sdm(is)=dm(i)
      else
       memb(i)=0
      endif
95    continue
      it=is
      iter=iter+1
      if(iter.lt.3) goto 91
c
c     computation of the parameters for all stars
c
      do 100 i=1,n
c      if(memb(i).eq.1) go to 100
        v0(i)=v(i)-3.2*ebmvm
        bmv0(i)=bmv(i)-ebmvm
        umb0(i)=umb(i)-eumbm 
        dm(i)=v0(i)-xmv0(i)
        adm(i)=dm(i)
100   continue
c
c     test for second iteration
c
      if(ipas.eq.0) then
       ipas=1
       goto 61
      endif
c
c     selection from the redenning
c
      do 101 i=1,n
      if(memb(i).eq.0) then
      ebdif(i)=ebmv2(i)-ebmvm
        if(ebdif(i).ge.0.10) memb(i)=8
c     print 800,i,ebmv2(i),ebdif(i)
c 800 format(2x,i3,2(2x,f5.2))
      endif
101   continue
c
c     opening of the output files on units 15 and 16
c
      inquire(file='f.out',exist=ll)
      if(.not.ll) then
      stat='new'
      else
      stat='old'
      endif
      open(15,file='f.out',status=stat)
c
      inquire(file='ubv.pem',exist=ll)
      if(.not.ll) then
      stat='new'
      else
      stat='old'
      endif
      open(17,file='ubv.pem',status=stat)
c
      inquire(file='evol.dat',exist=ll)
      if(.not.ll) then
      stat='new'
      else
      stat='old'
      endif
      open(16,file='evol.dat',status=stat)
c
      do 120 i=1,n
       write(16,510) i,sn(i),v0(i),dm(i),xmv0(i)
510    format(2x,i3,2x,i4,3(2x,f6.3))
120   continue
      write(6,520) ebmvm,i1,dmm
520   format(//,' e(b-v)mean = ',f6.3,' based on ',i3,' stars',//,
     1' m-m = ',f5.2,/)
105   format(2x,i3,1x,i4,2x,8(f6.3,2x),i1)
c
c     plot of the evolutionary difference curve
c
      call system("plt_evl")
      read(5,*) isu
      call system("ers")
      write(6,550)
550   format(/,' give the v0 limit:  ',$)
      read(5,*) v0l
      close(16,status='keep')
      open(16,file='evol.dat',status='old')
c
c     selection of members from the distance modulus
c
      ainf=-1.10
      asup=+1.55
      istar=0
      nn=n
      v0m=v0l-2.0
      do 1000 ij=1,9
      call moment(adm,nn,ddmm,adev4,sdev4,var4,skew4,curt4)
      do 108 i=1,n
      if(memb(i).gt.7) goto 108
      difdm(i)=ddmm-dm(i)
      if(difdm(i).lt.ainf) then
        memb(i)=9
        goto 108
        endif
      if(v0(i).lt.v0m) goto 108
      if(difdm(i).gt.asup) then
        memb(i)=9
        goto 108
      endif
      istar=istar+1
      adm(istar)=dm(i)
108   continue
      print *,ij,istar,' inf=',ainf,' sup=',asup,' ddmm=',ddmm
      ainf=ainf+0.1
      asup=asup-0.15
      nn=istar
      istar=0
1000  continue
      do 1001 i=1,n
      if(memb(i).lt.8) write(16,510) i,sn(i),v0(i),dm(i),xmv0(i)
1001  continue
c
c     computation of the evolutionary correction
c
      tebmv=0.0
      teumb=0.0
      tdm=0.0
      i1=0
      do 125 i=1,n
      if(v0(i).lt.v0l) then
      dif=dm(i)-dmg
c
      if(memb(i).eq.1.and.dif.lt.2.) then
       i1=i1+1
       tebmv=tebmv+ebmv(i)
       teumb=teumb+eumb(i)
       tdm=tdm+dm(i)
      endif
c
      if(dif.gt.3.5.and.umb0(i).lt.0.) then
      ebmve=0.0
      ebmv(i)=ebmve
      k=1
240   eumb(i)=0.70*ebmv(i)+0.05*ebmv(i)**2
      umb0(i)=umb(i)-eumb(i)
      do 258 li=1,14
      if(umb0(i).gt.ubsg(li).and.li.ne.1) goto 258
      if(li.eq.1) goto 257
      bmv0(i)=((bvsg(li)-bvsg(li-1))/(ubsg(li)-ubsg(li-1)))*
     1(umb0(i)-ubsg(li))+bvsg(li-1)
      goto 224
257   bmv0(i)=((bvsg(li)-bvsg(li+1))/(ubsg(li)-ubsg(li+1)))*
     1(umb0(i)-ubsg(li))+bvsg(li)
      goto 224
258   continue
224   if(k.gt.10) goto 124
      k=k+1
      ebmv(i)=bmv(i)-bmv0(i)
      v0(i)=v(i)-3.2*ebmv(i)
      goto 240
      endif
124   xmv0(i)=xmv0(i)+dif
      dm(i)=v0(i)-xmv0(i)
c     write(6,560) i,sn(i),v0(i),xmv0(i)
560   format(2x,i3,2x,i4,2(2x,f5.2))
      endif
c
      write(15,105) i,sn(i),v(i),bmv(i),umb(i),bmv0(i),umb0(i)
     1,v0(i),xmv0(i),dm(i),memb(i)
      if(memb(i).lt.8) then
      write(17,106) sn(i),v0(i),bmv0(i),umb0(i),xmv0(i),memb(i)
106   format(1x,i4,7x,f6.2,2x,f5.2,2x,f5.2,2x,f6.2,2x,i2)
c
      endif
125   continue
c
c     print of test value of the redenning without the most evolved stars
c
      ebmvm=tebmv/i1
      eumbm=teumb/i1
      dmm=tdm/i1
      write(6,850) ebmvm,eumbm,dmm,i1
850   format(/,' e(b-v) = ',f6.3,/,
     1' e(u-b) = ',f6.3,/,'    m-m = ',f5.2,/,
     2' with ',i3,' stars',/)
      close(16,status='keep')
      call system("plt_evl")
c
      stop
      end
c
c
      subroutine moment(data,n,ave,adev,sdev,var,skew,curt)
      dimension data(n)
      if(n.le.1)pause 'n must be at least 2'
      s=0.
      do 11 j=1,n
        s=s+data(j)
11    continue
      ave=s/n
      adev=0.
      var=0.
      skew=0.
      curt=0.
      do 12 j=1,n
        s=data(j)-ave
        adev=adev+abs(s)
        p=s*s
        var=var+p
        p=p*s
        skew=skew+p
        p=p*s
        curt=curt+p
12    continue
      adev=adev/n
      var=var/(n-1)
      sdev=sqrt(var)
      if(var.ne.0.)then
        skew=skew/(n*sdev**3)
        curt=curt/(n*var**2)-3.  
      else
        pause 'no skew or kurtosis when zero variance'
      endif
      return
      end