      program pktnor4
C     ==============================================
C        ver.   3-Mar-1996      S. Zola
C        popr.  13-April-2003
C        ver. 9-May-2003 also errors in flux added 
C        ver.15-Jan-2005 avoids flux errors =0.000 
C     ==============================================

      parameter (mdim=5000)
      implicit double precision (a-h,o-z)
      double precision magb,magv
      dimension phasB(mdim),magB(mdim),fazaB(mdim)
     1,phasV(mdim),magV(mdim),fazaV(mdim)
      character *10 fil1,fil2
      common /elem/xm0,p0
      common /param/x1,x2,x3,x4,s1,s2,fac
      common /f/factb,factv
      common /kolor/alfab,alfav,ss1,ss2
      common /range/ r1,r2,r3,r4
      common /odr/ px1,px2
      common /sig/czy,stala

C
C  data file 
C  Please give only up to 5 digits of JDhel or else this program will not
C  work  i.e. 50282.272727

      write(*,*) ' INPUT FILE: '
      read(*,'(a)') fil1
      write(*,*) ' OUTPUT FILE: '
      read(*,'(a)') fil2

      open(1,file=fil1)
      open(2,file=fil2)

      write(*,*) ' Threshold level for removing bad points?:  Nm,Nd '
C =============================================================================================================================
C          These parameters should be set by trial, when seen how many
C          single points is being removed from anormal point.
C  Nm - used while removing bad points near the minima
C  Nd - used while removing bad points outside the minima
C          They could be the same, though, if increase/decrease is fast,
C          I suggest they are to be different to account for fast brightness
C          change. Please note, removing bad pints is done only once!
C =============================================================================================================================
      read(*,*) px1,px2
C =============================================================================================================================
C  czy   -> factor1 = +-1. (To change sign of magnitudes if needed)
C  stala -> factor2  (To subtract a constant if delta_magnitudes are big)
C =============================================================================================================================
      read(1,*)czy,stala
      write(2,*)' factor =',czy,' subtracted: ', stala
C
C =============================================================================================================================
C  x1-x4 parameters describing phase ranges, s1,s2 - steps  in the ranges, where:
C  between  <0-x1> i <x4-1> normal points are calculated using step "s1" (primary minimum)
C  between  <x1-x2> i <x3-x4> normal points are calculated using step "s2"
C  between  <x2-x3> using step "s1*fac"  (secondary minimum)
C =============================================================================================================================
      read(1,*)x1,x2,x3,x4,s1,s2,fac
      if(x1.lt.-1.)stop
C =============================================================================================================================
C     coefficients for standard reduction
C     this is not done right, it is just a "zero-order" aproximation  
C     if the coefficients are set to 0 nothing is done (recommended)
C =============================================================================================================================
      read(1,*)alfab,alfav
      read(1,*)r1,r2,r3,r4
C  -------------------------------------------------------------------
C       Epoch (primary minimum)     Period
C  -------------------------------------------------------------------
      read(1,*)xm0,p0
C  -------------------------------------------------------------------
C  ===  factors to normalize LC in flux units for required phase      
C  The code must be run twice, set the coefficients to 1.0 for the initial
C  run and take the coefficients from the output file for the next run
C  -------------------------------------------------------------------
      read(1,*)factb,factv
      write(2,*) ' '
      write(2,*)'       M0           P '
      write(2,'(5x,f12.5,3x,f12.6)') xm0,p0
      write(2,*) ' '
      call czyt(phasb,magb,nobsb)
      call czyt(phasv,magv,nobsv)
      call pfaz(phasb,nobsb,fazab)
      call pfaz(phasv,nobsv,fazav)
      call sort3(fazab,phasb,magb,nobsb)
      call sort3(fazav,phasv,magv,nobsv)

      IF(ALFAb.ne.0..OR.ALFAV.NE.0.) THEN
      CALL systm(fazab,magb,nobsb,fazav,magv,nobsv)
      END IF

      write(2,*) '      ********************* '
      do 450 i=1,nobsb
      pr=(10.**(-magb(i)/2.5))/factb
450   write(2,460) i,phasb(i),fazab(i),magb(i),pr
      write(2,*) '      ********************* '
      do 451 i=1,nobsv
      pr=(10.**(-magv(i)/2.5))/factv
451   write(2,460) i,phasv(i),fazav(i),magv(i),pr
      write(2,*) '      ********************* '

      call sort(fazab,magb,nobsb,1)
      call sort(fazav,magv,nobsv,2)
      stop
460   format(1x,i4,2x,f12.5,3x,f8.5,5x,f6.3,3x,f6.4)
300   format(3x,f12.4,2x,f8.4,2x,f7.3)

      end

      subroutine czyt(phas,mag,nobs)
C
C Reading data - give only up to 5 digits of the JDhel, otherwise
C the program will not work
C
      parameter (mdim=5000)
      implicit double precision (a-h,o-z)
      double precision mag
      dimension phas(mdim),mag(mdim)
      common /sig/czy,stala
      do 74 i=1,5000
      read(1,*)phas(i),mag(i)
      mag(i)=(mag(i)*czy) + stala
      if(phas(i).lt.0.)stop
      if(phas(i)-60000.) 74,125,125
  74  continue
 125  nobs=i-1
      write(2,30) nobs
  30  format(/' liczba obserwacji =',i6/)
      return
      end

      subroutine pfaz(phas,nobs,faza)
      parameter (nbs=5000)
      implicit double precision (a-h,o-z)
C calculating phases with given M0 and P
      dimension phas(nbs),faza(nbs)
      common /elem/ xm0,p0
      do 100 i=1,nobs
      faza(i)=(phas(i)-xm0)/p0
      t=aint(faza(i))
      faza(i)=faza(i)-t
      if(faza(i).lt.0.) faza(i)=faza(i)+1.
 100  continue
      return
      end

      subroutine sort(faza,mag,nobs,ik)
      parameter (mdim=5000,ndm=5000)
      implicit double precision (a-h,o-z)
      double precision mag
      dimension xfazn(ndm),xmagn(ndm),xvar(ndm),numb(ndm)
      dimension xfaznfl(ndm),xfln(ndm),xvarfl(ndm),numbfl(ndm)
      dimension mag(nobs),faza(nobs),flux(nobs)
      dimension fl(mdim),weight(mdim)
      common /f/factb,factv
      factor=factb
      if(ik.eq.2)factor=factv
      call norm(faza,mag,nobs,xfazn,xmagn,numb,xvar,il)
      il=il-1
      write(2,55)
      do 501 i=1,il
      fl(i)=(10.**(-xmagn(i)/2.5))/factor
C     ----------------------------------------------
C      Next record assumes the weight for normal points for the 
C      DC W-D code  (here the same weight = 1. but can be set 
C      in a required way)
C     ----------------------------------------------
      weight(i)=1.
      write(2,54) xfazn(i),fl(i),xvar(i),xmagn(i),numb(i),i
  501 continue
      write(2,51)
      do 57 i=1,il,5
      write(2,56)xfazn(i),fl(i),weight(i),xfazn(i+1),fl(i+1),weight(i+1)
     $,xfazn(i+2),fl(i+2),weight(i+2),xfazn(i+3),fl(i+3),weight(i+3)
     $,xfazn(i+4),fl(i+4),weight(i+4)
57    continue
C =====================================================================
C     flux errors for testing only for now
C     howerver, it seems to work OK
C =====================================================================
      do 502 i=1,nobs
C      flux(i)=(10.**(-mag(i)/2.5)/factor)
       flux(i)=10.**(-mag(i)/2.5)
502   continue 
      call norm(faza,flux,nobs,xfaznfl,xfln,numbfl,xvarfl,iln)
      iln=iln-1
      write(2,*) '                                              '
      write(2,*) ' =============  FLUX errors =================='
      do 503 i=1,iln 
      if(xvarfl(i).lt.0.001) xvarfl(i)=0.001
      write(2,59) xfaznfl(i),xfln(i)/factor,xvarfl(i),numbfl(i),i
503   continue 
      write(2,*) ' =============  END FLUX errors =============='
      write(2,*) '                                              '
C =====================================================================
  52  format(i5)
  55  format(1x,'NORMAL POINTS (PHASE-FLUX-----MAGN.StDv.--NO of PTS'/)
  54  format(1x,f8.4,2x,f6.4,2x,f5.3,2x,f7.3,4x,i4,I5)
  51  format(1x,'N',4X)
  56  format(5(f6.4,1x,f6.4,f3.1))
  59  format(1x,f8.4,2x,f6.4,2x,f5.3,4x,i4,I5)
      RETURN
      END

      subroutine norm(faza,fmag,nobs,xfazn,xmagn,numb,xvar,iik)
      parameter (ndm=5000,nbs=5000)
      implicit double precision (a-h,o-z)
      dimension faza(nbs),fmag(nbs),xfazn(ndm),xmagn(ndm)
      dimension xfaz(ndm),xmag(ndm),numb(ndm),xvar(ndm)
      common /param/x1,x2,x3,x4,s1,s2,fac
      x0=0.
      iik=1
 888  continue
      if(x0.lt.x1.or.x0.ge.x4) then
      step=s1
      else
        if(x0.gt.x2.and.x0.lt.x3)then
        step=s1*fac
        else
        step=s2
        end if
      end if
      xk=x0+step
      if(xk.gt.1.05) go to 130
      k=0
      do 100 i=1,nobs
      if(faza(i).ge.x0.and.faza(i).lt.xk) go to 10
      go to 100
  10  k=k+1
      xfaz(k)=faza(i)
      xmag(k)=fmag(i)
 100  continue
      if(k.eq.0.) go to 999
      call sre(xfaz,k,xf)
      call sre(xmag,k,xm)
      if(k.eq.1) xvar(iik)=10.
      if(k.gt.1) call vari(xmag,k,var,xm)
      if(k.gt.2) call odrz(xfaz,xmag,k,var,xm,xf)
      if(k.gt.1) xvar(iik)=var
      if(k.eq.1) xvar(iik)=10.
      xfazn(iik)=xf
      xmagn(iik)=xm
      numb(iik)=k
      iik=iik+1
 999  x0=xk
      go to 888
130   return
      end

      subroutine vari(xmag,k,var,xm)
C  computing standard deviation
      parameter (ndm=5000)
      implicit double precision (a-h,o-z)
      dimension xmag(ndm),del(ndm)
      do 50 i=1,k
  50  del(i)=xmag(i)-xm
      d0=0.
      do 70 i=1,k
      d1=d0+del(i)*del(i)
      d0=d1
 70   continue
      var=d0/(float(k)*float(k-1))
      var=sqrt(var)
      return
      end

      subroutine systm(xb,yb,nb,xv,yv,nv)
      parameter (mdim=5000)
      implicit double precision (a-h,o-z)
C  ==========================================================
C  aproximate reduction to the standard system
C  WARNING!!!  IT may not work properly
C  alfa - coefficient
C  r1-r4 ranges             
C  ==========================================================
       dimension xb(nb),yb(nb),xv(nv),yv(nv),color(mdim)
       common /param/p1,p2,p3,p4,s1,s2,fac
       common /kolor/alfab,alfav,sss1,sss2
       common /rang/r1,r2,r3,r4
       iil=1
       x0=0.
999    continue
       if(x0.lt.r1.or.x0.ge.r4) then
       step=s1
       if(x0.gt.r2.and.x0.lt.r3) then
       step=fac*s1
       else
       step=s2
       end if
       end if
       xk=x0+step
       if(xk.gt.1.06) return
       call mean(xb,yb,nb,sr1,x0,xk,k1)
       call mean(xv,yv,nv,sr2,x0,xk,k2)
       color(iil)=sr1-sr2
       if(k1.eq.0.or.k2.eq.0) color(iil)=color(iil-1)
       do 101 i=1,nb
       if(xb(i).ge.x0.and.xb(i).lt.xk)yb(i)=yb(i)+color(iil)*alfab
  101  continue
       do 102 i=1,nv
       if(xv(i).ge.x0.and.xv(i).lt.xk)yv(i)=yv(i)+color(iil)*alfav
  102  continue
 104   x0=xk
       iil=iil+1
       go to 999

C    If there is a warning in your compiler about the next line, 
C    just ignore it       

       return
       end

      subroutine mean(p,fm,n,s,x0,xk,k)
      parameter (ndm=5000)
      implicit double precision (a-h,o-z)
      dimension p(n),fm(n),xmag(ndm)
      k=0
      do 100 i=1,n
      if(p(i).ge.x0.and.p(i).lt.xk) go to 10
      go to 100
   10 k=k+1
      xmag(k)=fm(i)
  100 continue
      if(k.eq.0) go to 350
      call sre(xmag,k,s)
      return
  350 s=30.
      return
      end

      subroutine sortr(x,y,n)
      parameter (ndm=5000)
      implicit double precision (a-h,o-z)
C   sorting
      dimension x(5000),y(5000)
c
      do 10 kk=1,n
      do 20 i=2,n
      xlo=dmin1(x(i-1),x(i))
      if(xlo.eq.x(i)) then
      tt=x(i-1)
      x(i-1)=xlo
      t=y(i-1)
      y(i-1)=y(i)
      x(i)=tt
      y(i)=t
      end if
  20  continue
 10   continue
      return
      end

      subroutine odrz(x,y,k,var,ysr,xf)
C removing bad points
      parameter (ndm=5000)
      implicit double precision (a-h,o-z)
      dimension x(k),y(k),x1(ndm),y1(ndm)
      common /param/p1,p2,p3,p4,s1,s2,fac
      common /odr/ px1,px2
      if(k.lt.2)return
      px=px1
      if(x(1).lt.p1.or.x(1).gt.p4) px=px2
      if(x(1).gt.p2.and.x(1).lt.p3) px=px2
      par=px*var
      l=1
      do 10 i=1,k
      p=abs(y(i)-ysr)
      if(p.gt.par) go to 9
      x1(l)=x(i)
      y1(l)=y(i)
      l=l+1
      go to 10
9     write(2,20)x(i),y(i),p,var,ysr
10    continue
      k1=l-1
      k=k1
      if(k1.eq.k) return
      if(k.gt.1) then
      call sre(x1,k1,xf)
      call sre(y1,k1,ysr)
      call vari(y1,k1,var,ysr)
      else
      xf=x1(1)
      ysr=y1(1)
      var=30.
      end if
20    format(1x,' pkt. wyrzucony: ',5(f8.4,1x))
      return
      end

      subroutine sre(x,n,sr)
      parameter (ndm=5000)
      implicit double precision (a-h,o-z)
C mean value od one dimensional matrix
      dimension x(ndm)
      a0=0.
      do 200 i=1,n
      a1=a0+x(i)
      a0=a1
200   continue
      sr=a0/float(n)
      return
      end

      subroutine sort3(x,y,z,n)
      parameter (ndm=5000)
      implicit double precision (a-h,o-z)
      dimension x(n),y(n),z(n)
c
      do 10 kk=1,n
      do 20 i=2,n
      xlo=dmin1(x(i-1),x(i))
       if(xlo.eq.x(i)) then
       tt=x(i-1)
       x(i-1)=xlo
       t=y(i-1)
       ttt=z(i-1)
       y(i-1)=y(i)
       z(i-1)=z(i)
       x(i)=tt
       y(i)=t
       z(i)=ttt
       end if
  20  continue
 10   continue
      return
      end
