c Program to fit a period-luminosity relation to Cepheid
c variables in the Large Magellanic Cloud and the Small
c Magellanic Cloud.  The basic data is in the two files
c lmc_cepheids.dat and smc_cepheids.dat.  You should
c use the V magnitudes in these files.  The contents of
c the columns in these files can be found at
c the top of columns.  Note that the formats are different
c in the two files.
c
c This program is incomplete.  It reads in and fits the LMC data
c only.  You need to read in and fit the SMC data, then derive
c the relative distance modulii of the SMC and LMC and compute
c the distance to the SMC assuming the LMC is at 51 kpc.
c
c Compile the program: f77 -o distance distance.f
c Run the program:     ./distance
c Examine the output:  more distance.out
c -------------------------------------------------------------
c
c Declare variables
c
      implicit double precision (a-h,o-z)
      real*8 vlmc(5000),plmc(5000),logplmc(5000)
c
c Open files
c
      open(unit=10,file='lmc_cepheids.dat',status='old')
      open(unit=50,file='distance.out')
      write(50,*) ' Assignment 3, Q.1'
      write(50,'(/,''Fits to the LMC and SMC Cepheid PL relations'')')
c
c Read in the LMC data
c
      do i=1,2
        read(10,'()')
      enddo
      do i=1,5000
        read(10,'(86x,f8.3,f12.7)',end=100) vlmc(i),plmc(i)
      enddo
 100  nlmc=i-1
c
c Convert P to log10P
c
      do i=1,nlmc
        logplmc(i)=log10(plmc(i))
      enddo
c
c Make a least squares fit to the LMC data V = almc*log10P + blmc
c
      call lsqfit(logplmc,vlmc,almc,blmc,nlmc)
      write(50,'(/,''LMC fit: V = '',f7.3,''*logP + '',f6.3)') almc,blmc
c
c Here the difference in distance modulus between the SMC and LMC is computed.
c This is done by evaluating the mean magnitudes of LMC and SMC Cepheids with 
c the approximately most-common period P=4 days and evaluating the difference.
c Also compute the SMC distance assuming the LMC is at 51 kpc.
c Write the results.
c
      v4lmc=almc*log10(4.0d0)+blmc
      write(50,'(''Mean V of LMC Cepheid with P=4d = '',f7.3)') v4lmc
c
c Finished
c
      stop
      end 
c
c
c
      subroutine lsqfit(x,y,a,b,n)
c Makes a least squares fit y=a*x+b to a set of n points (x_i,y_i)
c
c Declare variables
      implicit double precision (a-h,o-z)
      real*8 x(1),y(1)
c Compute some sums
      sx=0.0d0
      sy=0.0d0
      sxy=0.0d0
      sx2=0.0d0
      do i=1,n
        sx=sx+x(i)
        sy=sy+y(i)
        sxy=sxy+x(i)*y(i)
        sx2=sx2+x(i)*x(i)
      enddo
c Compute a and b
      d=sx*sx-n*sx2
      a=(sx*sy-n*sxy)/d
      b=(sx*sxy-sx2*sy)/d
c
      return
      end