c f77 aliens.f -o aliens
c This program computes the total mass of all aliens on Planet X, even
c though we can only detect aliens of mass greater than 30 kg.  Many
c simulations of this experiment are computed, and the standard
c deviation in the final result is our uncertainty.

      parameter (Nm=7,NmAll=10,dm=10.)
      parameter (Ntrials=50)
      real m(Nm),      N(Nm),      sig(Nm),     N_over_dm(Nm)
      real msim(NmAll),Nsim(NmAll),sigsim(NmAll)
      real mtotsim(Ntrials)
      data m/35,45,55,65,75,85,95/
      data N/66,54,46,34,26,14,5/
      data msim/5,15,25,35,45,55,65,75,85,95/
      idum=-479
      
      do j=1,Nm
       N_over_dm(j)=N(j)/dm
       sig(j)=sqrt(N_over_dm(j))
      enddo
      
      call fit(m,N_over_dm,Nm,sig,1,a,b,siga,sigb,chi2,q)
      print*,'slope and intercept:',b,a
      print*
      
      do k=1,Ntrials
       
       ! get inferred counts at 5,15,25, i.e., below our sensitivity
       do j=1,3
        Nsim(j)=b*msim(j)+a
        sigsim(j)=sqrt(Nsim(j))
        Nsim(j)=Nsim(j) + gasdev(idum)*sigsim(j)
       enddo
       ! now include (tweaks to) measured counts
       do j=4,NmAll
        Nsim(j)=N_over_dm(j-3) + gasdev(idum)*sig(j-3)
        sigsim(j)=sqrt(abs(Nsim(j))) ! in case < 0
       enddo
       do j=1,NmAll
        !print*,msim(j),Nsim(j),sigsim(j)
       enddo
       
       call fit(msim,Nsim,NmAll,sigsim,1,asim,bsim,siga,sigb,chi2,q)
       mtotsim(k)=0.
       do j=1,NmAll
        mtotsim(k)=mtotsim(k)+msim(j)*Nsim(j)*dm ! mult. by bin width
       enddo
       print 100,k,asim,bsim,mtotsim(k)
 100   format(i3,1x,f5.1,1x,f5.2,1x,f8.1)
       
      enddo
      
      print*
      call moment(mtotsim,Ntrials,ave,adev,sdev,var,skew,curt)
      print 110,'Total mass is',ave,'+/-',sdev
 110  format(a13,1x,f7.1,1x,a3,1x,f6.1)
      
c      ! compute total mass
c      total=0
c      total2=0
c      !do j=5,95,5
c      do j=0,100
c       total=total+(100.-j)*real(j)
c       total2=total2+real(j)*(100.-real(j))
c       !print*,j,total
c      enddo
c      print*,total,total2
      
      end ! of main program
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/ran1.f'
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/fit.f'
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/moment.f'
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/gasdev.f'
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/gammq.f'
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/gammln.f'
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/gser.f'
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/gcf.f'