*23456789a123456789b123456789c123456789d123456789e123456789f123456789g12
*
      program lc_fit
      parameter (npp=3,np=150,nxx=3,nx=1000)
      parameter (nyy=np*nxx,nyyy=nyy*nxx)
      dimension ip(npp),timep(np,npp),xmagp(np,npp),
     &          ix(nxx),timex(nx,nxx),xmagx(nx,nxx),
     &          xoffset(nxx),xmago(nxx),time(np)
      dimension iorder(nyy),iinteger(nyy),timexx(nyy),xmagxx(nyy,nxx)
      character cchar(nyy)*1
      character filenp(npp)*30,filenx(nxx)*30
      data filenp/'./lc_brunou.dat','./lc_brunob.dat',
     &            './lc_brunov.dat'/
*      data filenx/'./sn1972e/lc_u.dat','./sn1972e/lc_b.dat',
*     &            './sn1972e/lc_v.dat'/
*      data filenx/'./sn1984a/lc_u.dat','./sn1984a/lc_b.dat',
*     &            './sn1984a/lc_v.dat'/
*      data filenx/'./sn1986g/lc_u.dat','./sn1986g/lc_b.dat',
*     &            './sn1986g/lc_v.dat'/
*      data filenx/'./sn1981b/lc_u.dat','./sn1981b/lc_b.dat',
*     &            './sn1981b/lc_v.dat'/
*      data filenx/'./sn1990n/lc_u.dat','./sn1990n/lc_b.dat',
*     &            './sn1990n/lc_v.dat'/
*      data filenx/'./sn1992a/lc_u.dat','./sn1992a/lc_b.dat',
*     &            './sn1992a/lc_v.dat'/
      data filenx/'./sn1994d/lc_u.dat','./sn1994d/lc_b.dat',
     &            './sn1994d/lc_v.dat'/
*      data filenx/'./sn1994ae/lc_u.dat','./sn1994ae/lc_b.dat',
*     &            './sn1994ae/lc_v.dat'/
*      data filenx/'./sn1995d/lc_u.dat','./sn1995d/lc_b.dat',
*    &            './sn1995d/lc_v.dat'/
      data ix/nxx*0/,xmagxx/nyyy*0./
*
* The isuper is the year*1000 plus the order number of the letter
* designation.
*  a  b  c  d  e  f  g  h  i  j  k  l  m  n
*  1  2  3  4  5  6  7  8  9 10 11 12 13 14
*
*  SN 1972E isuper=1972000+5
*      isuper=1972000+5
*  SN 1981B isuper=1981000+2
*      isuper=1981000+2
*  SN 1984A isuper=1984000+1
*      isuper=1984000+1
*  SN 1986G isuper=1984000+1
*      isuper=1986000+7
* SN 1990N isuper=1990000+14
*      isuper=1990000+14
* SN 1992a isuper=1992000+1
*      isuper=1992000+1
* SN 1994d isuper=1994000+4
      isuper=1994000+4
* SN 1994ae isuper=1994000+31
*      isuper=1994000+31
* SN 1995d isuper=1995000+4
*      isuper=1995000+4
*
* Read in Bruno's template light curves.
*
      do 410 i=1,npp
      open(unit=1,file=filenp(i),status='old')
*
      call readhead(1,0,0)
*
      j=1
  120 continue
      read(1,*,end=130) timep(j,i),xmagp(j,i)
      if(isuper .eq. 1990000+14  .and.  i .eq. 3) then
* This is just because the V maximum of 90N occurred only
* one day later than the B maximum.  So we shift the V
* template to 1 day earlier.
          timep(j,i)=timep(j,i)-1
      end if
*      print*,j,timep(j,i),xmagp(j,i)
      j=j+1
      go to 120
  130 continue
      ip(i)=j-1
      close(unit=1)
*
  410 continue
*
*
* Read in the observed light curves.
*
      k=0
      do 420 i=1,nxx
      open(unit=2,file=filenx(i),status='unknown')
*
      j=1
      read(2,*,err=160,end=160)
      call readhead(2,0,0)
  150 continue
      read(2,*,err=160,end=160) timex(j,i),xmagx(j,i)
*      print*,timex(j,i),xmagx(j,i)
      if(isuper .eq. 1992000+1) then
          timex(j,i)=timex(j,i)-48622.7
        else if(isuper .eq. 1994000+4) then
          timex(j,i)=timex(j,i)-49432.90+18.
* This is Vacca & Leibundgut 1996 B max.
        else if(isuper .eq. 1994000+31) then
          timex(j,i)=timex(j,i)-9666.5
        else if(isuper .eq. 1995000+4) then
          timex(j,i)=timex(j,i)-9751.5
      end if
      k=k+1
      timexx(k)=timex(j,i)
      xmagxx(k,i)=xmagx(j,i)
*      print920,timex(j,i),xmagx(j,i)
*      print920,timex(j,i),0.,0.,xmagx(j,i),0.,0.
      j=j+1
      go to 150
  160 continue
      ix(i)=j-1
      close(unit=2)
*
  420 continue
*
      itimexx=k
*
      call ordering(2,itimexx,iorder,iinteger,timexx,cchar)
      do 421 i=1,nxx
      call ordering(-2,itimexx,iorder,iinteger,xmagxx(1,i),cchar)
  421 continue
      i=0
      do 422 k=1,itimexx
*      print920,k,timexx(k),(xmagxx(k,m),m=1,nxx)
      if(timexx(k) .gt. 0.) then
          i=i+1
          timexx(i)=timexx(k)
          do 423 m=1,nxx
          xmagxx(i,m)=xmagxx(k,m)
  423     continue
          do 424 l=k+1,itimexx
          if(timexx(l) .gt. timexx(i)+1.e-5) go to 170
          timexx(l)=-1.
          do 426 m=1,nxx
          if(abs(xmagxx(l,m)) .gt. 1.e-5  .and.
     &       abs(xmagxx(i,m)) .lt. 1.e-5       ) then
              xmagxx(i,m)=xmagxx(l,m)
            else if(abs(xmagxx(l,m)) .gt. 1.e-5) then
              xmagxx(i,m)=.5*(xmagxx(i,m)+xmagxx(l,m))
* I'm assuming there will never be any more than two
* separate observations with numerically identical 
* times.  If there are, the error is small:  I'm just
* doing an only slightly stupid unequal averaging.
          end if
  426     continue
  424     continue
  170     continue
      end if
  422 continue
      itimexx=i
*
      print*
      print908
  908 format('Photometry:  observed')
      do 428 j=1,itimexx
      print920,timexx(j),(xmagxx(j,i),i=1,nxx),0.,0.
  428 continue
      print930 
*
* Now find the mean offset between the templates and the
* observed light curves.
*
      print*
      do 430 i=1,nxx
      if(ix(i) .eq. 0) go to 430
      sum=0.
      icount=0
      do 440 j=1,ix(i)
*      if(isuper .eq. 1981000+2  .and.      ! This is for 
*     &  (timex(j,i) .lt. 14.    .or.       ! a maximum light fit. 
*     &   timex(j,i) .gt. 22.        )) go to 440
      if(isuper .eq. 1994000+4  .and.      ! This is for
     &  (timex(j,i) .lt. 40.    .or.       ! an early decline fit.
     &   timex(j,i) .gt. 100.        )) go to 440
      icount=icount+1
      call interpo(1,timex(j,i),timep(1,i),xmagp(1,i),
     &             ip(i),1,1,xmagpp)
*      print*,i,timex(j,i),xmagx(j,i),xmagpp
      sum=sum+(xmagx(j,i)-xmagpp)
  440 continue
      xoffset(i)=sum/real(icount)
      print*,'The ',i,'th offset magnitude is'
      print*,xoffset(i)
      print*
  430 continue
*
* Read in the full set of Bruno's times two late Branch & Doggett
* times.
*
      open(unit=3,file='./lc_brunot.dat',status='old')
      read(3,*) itime
      do 450 i=1,itime
      read(3,*) time(i)
  450 continue
*
* Print out the fitted light curve.
*
      print910
  910 format('Photometry:  fitted')
      do 460 i=1,itime
      do 470 j=1,nxx
      if(ix(j) .gt. 0) then
          call interpo(1,time(i),timep(1,j),xmagp(1,j),
     &          ip(j),1,1,xmagpp)
          xmago(j)=xmagpp+xoffset(j)
        else
          xmago(j)=0.
      end if
  470 continue
      print920,time(i),(xmago(j),j=1,nxx),0.,0.
*      print920,time(i),xmago(3)
  460 continue
*
  920 format(f10.3,5f10.4)
      print930
  930 format('END:')
*
      delx=1.
      do 480 i=1,nxx
      k=0
      do 482 j=1,itimexx
      if(abs(xmagxx(j,i)) .gt. 1.e-5) then
          k=k+1
          timex(k,i)=timexx(j)
          xmagx(k,i)=xmagxx(j,i)
          iorder(k)=j
      end if
  482 continue
      itimex=k
*
* smooth over 1 day.  Remember this means half DELX ahead
* and half DELX behind.
      call smooth(itimex,timex(1,i),xmagx(1,i),delx,0.)
      do 484 j=1,itimex
      xmagxx(iorder(j),i)=xmagx(j,i)
      print920,timex(j,i),xmagx(j,i)
  484 continue
*
  480 continue
*
      print*
      print940,delx
  940 format('Photometry:  observed:  box-car smoothed',
     & ' over a ',f5.2,' day interval')
      do 500 j=1,itimexx,1
      print920,timexx(j),(xmagxx(j,i),i=1,nxx),0.,0.
  500 continue
      print930

*
      end
