C ====================================================================
C     Version: 6-Sep-1988    G.P.
C             13-Sep-1988    new version
C              2-Dec-1993    Lahey Fortran version (S.Z.)
C             22-Oct-1994    some bugs corrected
C             25-Mar-1996    I filter added
C              8-gru-2005    poprawka starych algorymtmow  (S.Z)
C ====================================================================
      PROGRAM FinalRed
      implicit double precision(a-h,o-z)
      character*60 line0,line1,line2
      character*30 FInp,FOut
      character*30 SName
      character *4 filtry(6)
      character*3 Month,Mnth1
      character*1 key,ft1(10)
      character *1 filter
      parameter (ndim=5000)
      dimension time(ndim),flux(ndim),dmag(ndim),phase(ndim),dme(ndim)
      dimension zv(ndim)
c  *******************************
c     dimension ndni(12),nmr(12)
c     character*3 mies(12)
c     data mies/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP',
c    $'OCT','NOV','DEC'/
c     data /nmr/1,2,3,4,5,6,7,8,9,10,11,12/
c     data ndni/31,28,31,30,31,30,31,31,30,31,30,31/
c     equivalence (mies(1),nmr(1)),(mies(2),nmr(2)),(mies(3),nmr(3)),
c    $(mies(4),nmr(4)),(mies(5),nmr(5)),(mies(6),nmr(6)),
c    $(mies(7),nmr(7)),(mies(8),nmr(8)),(mies(9),nmr(9)),
c    $(mies(10),nmr(10)),(mies(11),nmr(11)),(mies(12),nmr(12))
c  *******************************
      write(*,*)' INPUT Data File Name: '
      read(*,'(a)')FInp
      write(*,*)' OUTPUT Data File Name: '
      read(*,'(a)')FOut
      filtry(1)='U'
      filtry(2)='B'
      filtry(3)='V'
      filtry(4)='R'
      filtry(5)='I'
      filtry(6)='-'
C
C***************** Input Data of MAIN & COMP Stars *********************
C
      call TakeCoor(SName,xm0,p,av,dv,ac,dc)
      XmType=xm0+2400000.
      call RAConv(av)
      call DCConv(dv)
      call RAConv(ac)
      call DCConv(dc)
      open(1,file=FInp)
      open(2,file=FOut)
      open(3,file='tst1.fredm')
C     open(9,file='tst2.fredn')
C
C******************** General DATA Input *******************************
C
      write(*,*)'Mean Extinction Coefficients (Y/N) [Y]: '
      read(*,'(a)')KEY
      write(*,'(1(/1x))')
      read(1,'(a)')line0
77    ntot=1
      MODE=1
      mft=0
      read(1,'(a)',END=999,ERR=1000)line1
      write(*,'(1x,a60)') line1
      read(1,763)line2
763   format(A60)
      write(*,'(1x,a60)') line2
      read(line1,'(6x,i2,1x,a3,1x,i4,32x,i2)')NDay,Month,NYear,NZone

C  ================    To jest tylko dla testow !!! ===========

      ihg0=50
      do 345 ihg=1,10
      ft1(ihg)=line2(ihg0+ihg:ihg0+ihg)
345   continue

      do 10 iip=1,10
10    if(ft1(iip).ne.' ') filter=ft1(iip)
      if(ft1(7).ne.' ') then
      filter=ft1(7)
      else
      filter=ft1(1)
      end if
C      write(*,'(10x,a1)') filter

C  ================    To ^ jest tylko dla testow !!! ===========

      do 654 iut=1,6
      if(filter.eq.filtry(iut)) then
      mft=1
      end if
654   continue

      if(mft.eq.0) then
      write(*,*) ' This filter name N/A - Enter new filter name: '
      read(*,'(a1)') filter
      end if

      call ExtCf(filter,aext,key)
      write(*,765) filter
765   format(5x,25HReducing data in filter: ,A1)
33    read(1,*,ERR=555,END=555)ndat,time(ntot),flux(ntot)
      ntot=ntot+1
      goto 33
C
C************** END of MAIN DATA READING *******************************
555   ntot=ntot-1
      write(*,*) '     Number of data: ',ntot
C************** START of REDUCTION *************************************
C
      NMon=MConv(Month)
      call JD(NDay,NMon,NYear,hjd,tit)
C     write(9,'(" tint : ",F12.5)') tint
C     write(9,'(" hjd  : ",F12.5)') hjd

      if(MODE.eq.1) then
      BACK=-1
      if(time(1).gt.11.) BACK=-BACK
      HJD1=hjd+BACK
      call JDConv(HJD1,ND1,Mnth1,NY1)
      endif

      STEP=(time(ntot)+time(1))/48.
      tt=tit+STEP/36525.
c     call HelCorr(av,dv,tt,hc)
      write(3,*) 'HJD1   before   call  hellcorrn, HJD,tt ', HJD1,HJD,tt
      call HelCorrn(av,dv,hjd1,hc)
C     write(9,*) 'hjd1', HJD1
      do 11 i=1,NTOT
      time(i)=(time(i)-NZone)/24.
      if(flux(i).gt.0.) then
      dmag(i)=-2.5*dlog10(flux(i))
      else
      dmag(i)=77.777
      end if
11    continue
C     tgr=SidGr0(tit)
      tgr=SidGr0n(HJD,tt)
      write(3,'(" tgr: ",F14.8)')  tgr
      call MainExt(time,av,dv,ac,dc,dmag,dme,tgr,aext,zv,NTOT)
      do 200 k=1,ntot
C     write(9,'(" time:",F12.5," HJD: ",F12.5, " hc: ", F8.5)') 
C    $ time(k),hjd,hc
      time(k)=time(k)+hjd+hc
C     write(9,'(" time:",F12.5," HJD: ",F12.5, " hc: ", F8.5)') 
C    $ time(k),hjd,hc
C     write(9,*) 'hjd: ', hjd
C     write(9,*) 'hc: ',  hc

200   continue
      call CPhase(xm0,p,time,phase,NTOT)
C
C************** END of One Filter REDUCTION ****************************
C*************** Write Data to OUTPUT File *****************************
C
C     if(MODE.gt.1) goto 111

c  ----------------------------------------------------------------------
c  ----------------------------------------------------------------------
      if(BACK.eq.1.) then
c  ----------------------------------------------------------------------
      write(2,'(1x,''DATE: '',i2,''-'',a3,''-'',i4,''/'',7x,''OBJECT: ''
     &,a30)')Nday,Month,NYear,SName
        IF(MONTH.EQ.'FEB') THEN
        nday1=NDAY+1
        if(Nday.eq.28) nday=1
        ELSE
        IF(MONTH.EQ.'APR'.OR.MONTH.EQ.'JUN'.OR.MONTH.EQ.'SEP'.OR.
     &   MONTH.EQ. 'NOV') THEN
        nday1=NDAY+1
        if(Nday.eq.30) nday1=1
        ELSE
        nday1=NDAY+1
         if(Nday.eq.31) nday1=1
         END IF
         END IF

      write(2,'(6x,''/'',i2,''-'',a3,''-'',i4,2x,''JD0= '',f12.4,2x,
     &''P= '',f17.10,''  k='',f5.2)')NDay1,Mnth1,NY1,XmType,p,aext

      else

      write(2,'(1x,''DATE: '',i2,''-'',a3,''-'',i4,''/'',7x,''OBJECT: ''
     &,a30)')ND1,Mnth1,NY1,SName
      write(2,'(6x,''/'',i2,''-'',a3,''-'',i4,2x,''JD0= '',f12.4,2x,
     &''P= '',f17.10,''  k='',f5.2)')Nday,Month,NYear,XmType,p,aext
c  ----------------------------------------------------------------------
      endif
c  ----------------------------------------------------------------------
c  ----------------------------------------------------------------------

      write(2,'(1x,67(''-''))')
      write(2,'(3x,''JDhel'',8x,''DMag'',5x,''Phase'',4x,
     &''ExtCorr'',2x,''ZDist'',4x,''FILTER'',2x,''No.'')')
      write(2,'(1x,5(''-''),'' 2400000.+ '',51(''-''))')
111   MODE=MODE+1
      do 88 i=1,NTOT
      write(2,'(1x,f11.5,2x,f7.3,1x,f9.5,3x,f6.3,3x,f4.1,7x,a1,2x,I5)')
     &time(i),dme(i),phase(i),dmag(i)-dme(i),zv(i),filter,i
88    continue
      write(2,'(1x,67(''-''))')
C
C*************** Next FILTER Reduction *********************************
C
      backspace(1)
      goto 77
1000  write(*,*)'  ERROR OCCURED READING INPUT FILE !!! '
      goto 1001
999   continue
      close(2)
      close(1)
1001  STOP
      end

      FUNCTION MConv(Month)
      character*3 Month,mon(12)
      data mon/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP',
     &'OCT','NOV','DEC'/
      do 1 i=1,12
      if(Month.eq.mon(i)) then
      mconv=i
      return
      endif
1     continue
      return
      end

      SUBROUTINE HelCorrn(Alf,Dec,hjd,HC)
      implicit double precision(a-h,o-z)
C
C  Sun coordinates calculated according to aproximate forlmulae from the
C      Astronomical Almanac  (S. Zola,  May 2005)
C
      pi=4.d0*DATAN(1.d0)
      xconv = pi/180.d0
      Xn = HJD - 51545.d0
C     Xn = Xn-1
      write(3,*)' in helcor, before call SunCoorn,  HJD, Xn: ',  HJD, Xn
      eps = 23.439 -0.0000004 * Xn
      factor=0.0057755
      eps = eps * xconv
      call SunCoorn(HJD,Slon,RSun)
C     write(9,*) " Slon, eps, Rsun ", slon,eps,rsun
c      hc=dcos(alf)*dcos(dec)*dcos(slon)
c      h0=(dsin(alf)*dcos(dec)*dcos(eps)+dsin(dec)*dsin(eps))*dsin(slon)
c      hc=hc+h0
c      hc=-hc*rsun*factor
       hc1 = dcos(slon)*dcos(alf)*dcos(dec)
       hc2 = dsin(eps)*dsin(dec)+dcos(eps)*dsin(alf)*dcos(dec)
       hc2 = hc2 * dsin(slon)
       hc = factor*Rsun*(Hc1+hc2)
       hc = -hc
      write(3,*) '  ======= In  HelCorr ========='
      write(3,*) '  hjd        Lambda_Slon         Rsun  '
      write(3,'(1x,F12.6,1x,F12.7,2x,F10.7)') HJD,Slon,Rsun
      write(3,'(14x,F12.7,2x,F10.7)') Slon/pi*180.,Rsun
      write(3,*) '    hel. corr.      eps      '
      write(3,'(1x,F12.6,2x,F10.7)')  hc,eps
      write(3,*) '    RA             Dec     '
      write(3,'(1x,F12.6,1x,F12.6)')  alf,dec
      write(3,'(1x,F12.6,1x,F12.6)')alf/pi*180.d0/360.*24.,dec/pi*180.d0
      return
      end

      SUBROUTINE RAConv(alfa)
C*************** Change HH.MMSS-->RAD **********************************
      implicit double precision(a-h,o-z)
      pi=4.*datan(1.d0)/180.
      h=dint(alfa)
      h0=(alfa-h)*100.
      hm=dint(h0)
      hs=(h0-hm)*100.
      alfa=((hs/60.+hm)/60.+h)*15.*pi
      return
      end

      SUBROUTINE DCConv(dec)
C*************** Change oo.''""-->RAD **********************************
      implicit double precision(a-h,o-z)
      pi=4.*datan(1.d0)/180.
      SGN=1.
      if(dec.lt.0.) SGN=-1
      h=dint(SGN*dec)
      h0=(SGN*dec-h)*100.
      hm=dint(h0)
      hs=(h0-hm)*100.
      dec=((hs/60.+hm)/60.+h)*pi*SGN
      return
      end


      SUBROUTINE JD(NDay,NMonth,NYear,HJD,TInt)
      implicit double precision(a-h,o-z)

      if(NMonth.le.2) then
       qmon=NMonth+12.
       year=NYear-1.
      else
       qmon=NMonth
       year=NYear
      endif
      a=dint(year/100.)
      aa=dint(a/4.)
      b=2.-a+aa
      bb=dint(365.25*year)+dint((30.6001*(qmon+1)))
      HJD=bb+b+NDay-679005.5
      TInt=(hjd-15020.)/36525.
      return
      end

      SUBROUTINE SunCoorn(HJD,SLon,RSun)
C
C   introduced in May 2005 -  aproximate formulae from the
C       Astronomical Almanac     (S. Zola)
C
      implicit double precision (A-H,O-Z)
      pi=4.d0*DATAN(1.d0)
      pi2= 2.d0* Pi
      xconv = pi/180.d0
      Xn = HJD - 51545.d0
C     Xn = Xn -1.
      eps = 23.439 -0.0000004 * Xn
      Gdeg = 357.528 + 0.9856003 * Xn
      XL = 280.460 + 0.9856474 * Xn
      write(3,*)  '  ------------- in  SunCoorn  ---------------------'
      do 133 ikl=1,1000
      if(XL.lt.0.) Xl = Xl + 360.d0
      if(Xl.gt.360.d0) Xl = Xl - 360.d0
      if(Gdeg.lt.0.) Gdeg = Gdeg + 360.d0
      if(Gdeg.gt.360.d0) Gdeg = Gdeg - 360.d0
133   continue
      write(3,*)  '  Gdeg           XL '
      write(3,'(1x,F12.6,1x,F12.6)')  Gdeg,Xl
      G = Gdeg * xconv

      Slon = Xl + 1.915d0*dsin(G) + 0.020d0*dsin(2.d0*G)
      write(3,*)  '  Slon      '
      write(3,'(1x,F12.6)')  slon
C     if(Slon.lt.0.d0) Slon = Slon + 360.
C     if(Slon.gt.360.d0) Slon = Slon - 360.
      Slon = Slon * xconv
      RSun= 1.00014 - 0.01671 * dcos(G) - 0.00014*dcos(2.d0*G)
      write(3,*)  '  -------------------------------------------------'
      return
      end

      SUBROUTINE SunCoor(T,SLon,RSun)
      implicit double precision (A-H,O-Z)
C  proceudra juz zle liczy wspolrzedne Slonca
C
C     pii=3.141592653d+00
      pii=4.*datan(1.d0)
      pi=pii/180.d0
      ch=2.*pii
      ql=2.8122183333d2+1.719174999d0*t+4.527777d-4*t**2+3.33d-6*t**3
      qm=3.584758333d2+3.599904975277d4*t-0.00015*t**2-3.333d-6*t**3
      ql=ql*pi
      qm=qm*pi
      qm=dmod(qm,ch)
      ql=dmod(ql,ch)
      em=.01675104-4.18d-5*t-1.26d-7*t**2
      cc=(2.*em-.25*em**3+(5./96.)*em**5)*dsin(qm)
      cc=cc+(1.25*em**2-(11./24.)*em**4)*dsin(2.*qm)
      cc=cc+((13./12.)*em**3-(43./64.)*em**5)*dsin(3.*qm)
      cc=cc+((103./96.)*em**4)*dsin(4.*qm)
      v=qm+cc
      xl=v+ql
      ro=((1.-em**2)/(1.+em*dcos(v)))*1.00000023
      a=pi*(153.23+2.25187541d4*t)
      b=pi*(216.57+4.50375082d4*t)
      c=pi*(312.69+3.29643577d4*t)
      d=pi*(350.74+4.452671142d5*t-.00144*t**2)
      e=pi*(231.19+20.20*t)
      h=pi*(353.4+6.59287155d4*t)
      a=dmod(a,ch)
      b=dmod(b,ch)
      c=dmod(c,ch)
      d=dmod(d,ch)
      e=dmod(e,ch)
      h=dmod(h,ch)
      xx=pi*(.00134*dcos(a)+.00154*dcos(b)+.002*dcos(c))
      xx=xx+pi*(.00179*dsin(d)+.00178*dsin(e))
      xy=5.43d-6*dsin(a)+1.575d-5*dsin(b)+1.627d-5*dsin(c)
      xy=xy+3.076d-5*dcos(d)+9.27d-6*dsin(h)
      xle=xl+xx
      SLon=dmod(xle,ch)
      RSun=ro+xy
      return
      end

      SUBROUTINE ExtCf(Filter,aext,key)
      implicit double precision(a-h,o-z)
      character*1 Filter,f(5),KEY
      dimension COEF(5)
      data f/'U','B','V','R','I'/
c  --------------------------------------------------------
c  dopisalem wsp. ekst. dla filtrow R & I (bez wyznaczania)
c  --------------------------------------------------------
      data COEF/0.77,0.40,0.28,0.25,0.20/
      if(KEY.eq.'N'.or.KEY.eq.'n') then
      write(*,'(a)') 'Ext. Coeff. for Filter: ',Filter
      read(*,*) aext
      return
      endif
      do 3 i=1,5
      if(Filter.eq.f(i)) then
      aext=COEF(i)
      return
      endif
3     continue
      aext=COEF(2)
      return
      end

      SUBROUTINE CPhase(xm0,Per,hjd,phase,NTOT)
      implicit double precision(a-h,o-z)
      dimension hjd(ntot),phase(ntot)
      do 10 i=1,NTOT
      phase(i)=(hjd(i)-xm0)/Per
      phase(i)=phase(i)-dint(phase(i))
      if(phase(i).lt.0.) phase(i)=phase(i)+1.
10    continue
      return
      end

      SUBROUTINE MainExt(time,av,dv,ac,dc,dmag,dme,tgr,aext,zv,NTOT)
      implicit double precision(a-h,o-z)
      dimension time(NTOT),dmag(NTOT),dme(NTOT),zv(NTOT)
C
C            Lattitude and longitude should be in radians
C
C     data SFi,SLag/0.872624557, -0.374877337/                          KRK50cm
      data SFi,SLag/0.86539242,-0.35066573/                             Suhora
C     data SFi,SLag/0.66275966,-0.39482266/                             Kryonerion
      pi=4.*datan(1.d0)/180.
      Q= 1.00273790935
      CONS=0.0012
      cons1=0.0018167
      cons2=0.002875
      cons3=0.0008083
      do 100 i=1,NTOT
      write(3,*) ' TGR, time in MainExt: ',tgr*24.,time(i)*24.
      tv=(tgr+Q*time(i))*24.*15.*pi-SLag
      tgv=tv-av
      tgc=tv-ac
      szv=1./(dsin(SFi)*dsin(dv)+dcos(SFi)*dcos(dv)*dcos(tgv))
      szc=1./(dsin(SFi)*dsin(dc)+dcos(SFi)*dcos(dc)*dcos(tgc))
      cosv=1./szv
      zv(i)=dacos(cosv)/pi
C      secz â^À^Ó 0.0018167(sec z -1) â^À^Ó  0.002875(sec z -1)Â? â^À^Ó 0.0008083(sec z â^À^Ó1)Â? 
C    stare przyblizenie na airmass
C     xv=szv*(1.-CONS*(szv*szv-1.))
C     xc=szc*(1.-CONS*(szc*szc-1.))
C
C    New aproximation for airmass
C
      xv=szv-cons1*(szv-1.)-cons2*(szv -1.)**2-cons3*(szv-1.)**3
      xc=szc-cons1*(szc-1.)-cons2*(szc -1.)**2-cons3*(szc-1.)**3
      dx=xv-xc
      if(xv.lt.4.0)then
      de=dx*aext
      else
      de=999999.
      endif
      dme(i)=dmag(i)-de
100   continue
      return
      end


      FUNCTION SidGr0n(HJD0,TT)
      implicit double precision(a-h,o-z)
C
C  Procedure introduced in May 2005 after prescription taken
C  from  http://aa.usno.navy.mil/faq/docs/GAST.html
C
      D = HJD0 + TT - 51545.0
      D0 = HJD0 - 51545.0
      T = D/3652
      GMST = 6.697374558 + 0.06570982441908* D0 +  0.000026 * T*T
       do 10 i = 1,20
       if(GMST.GT.24.d0) GMST = GMST - 24.
       if(GMST.LT.0.d0) GMST = GMST + 24.
10     continue
       SidGr0n = GMST /24.
      return
      end

      FUNCTION SidGr0(TInt)
      implicit double precision(a-h,o-z)
C  old folrmula by G. Pajdosz
      T=6.6460655+2400.051262*TInt+2.58622d-5*TInt*TInt
      if(T.lt.0.0) T=T+2400.0
      SidGr0=(T-24.0*dint(T/24.0))/24.
      return
      end

      SUBROUTINE JDConv(hjd,NDay,Month,NYear)
      implicit double precision (a-h,o-z)
      character*3 Month,mon(12)
      data mon/'JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP',
     &'OCT','NOV','DEC'/
      z=dint(hjd+2400000.5)
      a=(dint(z-1867216.25))/36524.25
      aa=dint(a/4.)
      b=z+a-aa+1525.
      c=dint((b-122.1)/365.25)
      d=dint(365.25*c)
      f=dint((b-d)/30.6001)
      ee=dint(30.6001*f)
      Nday=b-d-ee
      if(f.lt.13.5) qmon=f-1.
      if(f.ge.13.5) qmon=f-13.
      if(qmon.lt.2.5) NYear=c-4715.
      if(qmon.ge.2.5) NYear=c-4716.
      NMon=qmon
      Month=mon(NMon)
      if(Nday.eq.0) then
      hjd=hjd-1.
      z=dint(hjd+2400000.5)
      a=(dint(z-1867216.25))/36524.25
      aa=dint(a/4.)
      b=z+a-aa+1525.
      c=dint((b-122.1)/365.25)
      d=dint(365.25*c)
      f=dint((b-d)/30.6001)
      ee=dint(30.6001*f)
      Nday=b-d-ee+1.
      if(f.lt.13.5) qmon=f-1.
      if(f.ge.13.5) qmon=f-13.
      if(qmon.lt.2.5) NYear=c-4715.
      if(qmon.ge.2.5) NYear=c-4716.
      NMon=qmon
      Month=mon(NMon)
      endif
      return
      end

      SUBROUTINE TakeCoor(Name,xm0,Per,av,dv,ac,dc)
      implicit double precision(a-h,o-z)
      parameter (nst=500)
      dimension am(nst),p(nst),a1(nst),d1(nst),a2(nst),d2(nst)
      character*30 OldNam(nst),Name
      character*73 LINE
      character*1 key
C  ================================================================
C     a full path to the catalogue with data of stars declared here
C  ================================================================
      open(7,file='./fcatal.str')
19    rewind 7
      IFAT=0
      write(*,*)' Enter STAR NAME [CK Boo]: '
      read(*,'(a)')Name
      do 9 l=1,nst
      read(7,'(a)',END=500)OldNam(L)
      read(7,'(a)')LINE
      backspace(7)
      read(7,*)am(l),p(l),a1(l),d1(l),a2(l),d2(l)
      if(Name.eq.OldNam(l)) then
      xm0=am(l)-2400000.
      Per=p(l)
      av=a1(l)
      dv=d1(l)
      ac=a2(l)
      dc=d2(l)
      write(*,*)'     JD0              Period     R.A.var     Dec.var
     & R.A.cmp    Dec.cmp'
      write(*,'(a)')LINE
      goto 600
      endif
9     continue
600   write(*,*)'Edit, Choose Another, Accept (E/C/A) [A]: '
      read(*,'(a)')key
      if(key.eq.'c'.or.key.eq.'C') goto 19
      if(key.eq.'e'.or.key.eq.'E') then
      IFAT=1
      goto 700
      endif
      close(7)
      return
500   continue
      OldNam(l)=Name
      write(*,*)'   NO SUCH STAR !!!'
700   write(*,*)'   Enter JD0 , PERIOD: '
      read(*,*)xm0,Per
      am(l)=xm0
      xm0=xm0-2400000.
      p(l)=Per
      write(*,*)' MAIN STAR - R.A.(hh.mmss): '
      read(*,*)av
      a1(l)=av
      write(*,*)'             Dec.(oo.``""): '
      read(*,*)dv
      d1(l)=dv
      write(*,*)' COMP STAR - R.A.(hh.mmss): '
      read(*,*)ac
      a2(l)=ac
      write(*,*)'             Dec.(oo.``""): '
      read(*,*)dc
      d2(l)=dc
      write(*,*)'   SAVE THIS DATA (Y/N) [Y]: '
      read(*,'(a)')key
      if(key.eq.'n'.or.key.eq.'N') goto 19
      if(IFAT.eq.1)then
      do 29 k=l+1,nst
      read(7,'(a)',END=555)OldNam(k)
      read(7,*)am(k),p(k),a1(k),d1(k),a2(k),d2(k)
29    continue
555   l=k-1
      endif
      rewind 7
      do 99 i=1,l
      write(7,'(a)')OldNam(i)
      write(7,'(1x,f13.5,2x,f15.10,2x,f8.5,2x,f9.5,2x,f8.5,2x,f9.5)')
     &am(i),p(i),a1(i),d1(i),a2(i),d2(i)
99    continue
      close(7)
      return
      end