#!/usr/local/bin/perl
#(C) Joseph Mack NA3T 2000, jmack@wm7d.net 
#code released under GPL.
#Code from article on the Aurora opening of 15-16 Jul 2000, 
#presented at the N.E.W.S conference 26-28 Aug 2000, Enfield CT.
#
use strict;
$|=1;

#simulate_4.pl:
#help given if incorrect number of parameters given
#monte carlo added
#calculation of variances fixed up.
# 
#simulate_3.pl is a modified version of simulate_2.pl
#simulate_2.pl only ran 1 simulation and produced an error estimate for the fit
#however since the simulation is probablistic, it is hard to tell when the fit
#is getting close to convergeance.
#simulate_3.pl runs several fits, and produces an average fit and variances.

#how many times should you expect to see each call of a population of n operators, 
#if you make a certain number of contacts

#operators do not have the same probability of being worked.
#the probability is set by a function or by manually set bins.
#this is to make the probability density function an exponential
#rather than the concave downward of the standard poisson.

=head1 NAME 

simulate_4.pl - determines fit for a given population size and a scaling_factor to 
the contact frequency distribution of a ham radio contest or propagation opening.
The derivation of the model and the use of this program are described in the report
on the Aurora opening of 15-16 Jul 2000 on http://www.wm7d.net/azproj.shtml.

=head1 AUTHOR

Joseph Mack (Joe NA3T), jmack@wm7d.net

=head1 COPYRIGHT

(C) Joseph Mack 2000

=head1 DISTRIBUTION/LICENSE

GPL

=head1 SYNOPSIS

B<./simulate_4.pl> population scaling_factor

eg

./simulate_4.pl 800 0.025

The population is the expected number of people in the opening. The scaling_factor
is a number describing the distribution of probabilities of making a contact for the
operators in the opening. The program will output the statistics for the fit. The
best fit is determined from multiple runs of this program with different arguements.

=head1 HARDCODED INFORMATION

The code contains the model for the distribution of probabilities of each ham making a contact.
This is in the routine offer_contacts. The only code left is case 2, the other models having been
excluded.

The code contains the number of contacts that were seen in the logs

my $contacts = 2598;

and the frequence distribution seen in the logs

my @data= ("0","184", "113", "63", "48", "34", "18", "16", "20", "12", "8", "10", "5", "7", "4", "4",
"6", "2", "5", "1", "1", "2", "1", "1", "2", "1", "1", "1", "1", "1", "1" );

Alternate versions of @data can be entered for monte carlo estimation of errors.

The program produces a probablistic fit. 
Usually large numbers (64-1024) of runs are done to find the mean and standard
deviation of the best fit. The number of runs is set by 

my $simulations=1024

=cut
#------------------------------

#Uncomment one of these sets
#1. real contest set
#2598 contacts, 571 calls
#the number of people who were worked 0 times is unknown (here entered as 0).
my $contacts = 2598;
#experimental data
my @data= ("0","184", "113", "63", "48", "34", "18", "16", "20", "12", "8", "10", "5", "7", "4", "4", "6", "2", "5", "1", "1", "2", "1", "1", "2", "1", "1", "1", "1", "1", "1" );
#monte carlo data
#00
my @data=("0","183.1","114.0","66.5","56.2","34.9","17.1","15.4","17.1","7.7","5.1","11.1","5.5","6.8","4.0","3.4","6.6","0.9","4.0","2.3","0","2.6","0","2.3","2.7","0.9","1.6");
#01
#my @data=("0","168.2","107.6","69.2","48.7","36.7","25.3","21.3","20.9","11.4","3.6","5.9","5.6","11.7","2.2","2.6","6.2","1.8","2.5","1.0","0","3.1","0.7","0.8","2.4","3.1","0.1");
#02
#my @data=("0","187.9","120.7","60.1","56.0","35.7","14.6","9.4","19.8","9.8","11.1","10.8","7.5","7.1","0.5","0","6.0","0.3","5.3","3.6","0","3.8","0.9","1.5","2.1","1.3","2.5");
#03
#my @data=("0","178.9","110.8","62.6","42.4","28.5","15.5","15.7","18.6","20.1","6.7","5.8","0","7.2","4.7","5.6","5.1","0.7","1.3","2.1","0.8","0.6","2.0","1.6","1.2","2.0","1.0");
#04
#my @data=("0","204.0","117.8","63.3","57.3","33.2","20.1","14.4","20.8","12.1","9.3","10.9","1.3","6.2","6.5","4.5","4.2","4.6","7.5","0.9","0.6","0.7","1.5","2.7","2.9","0.6","1.8");
#05
#my @data=("0","188.8","97.2","61.6","51.3","34.1","18.4","12.5","30.9","8.7","6.6","6.4","4.9","10.5","6.6","5.4","8.7","1.6","4.3","0","0.3","0.5","2.4","0.7","1.5","1.7","2.3");
#06
#my @data=("0","181.8","101.1","62.0","51.9","27.9","17.5","13.9","20.6","10.5","5.3","10.6","5.1","3.2","4.0","2.4","4.8","2.7","6.6","3.8","0","2.0","1.6","1.6","1.4","1.0","1.2");
#07
#my @data=("0","178.8","117.5","40.4","36.4","29.1","17.6","17.8","18.8","11.0","7.8","10.9","2.0","6.0","2.5","2.8","6.0","4.9","2.4","1.8","2.5","3.5","1.2","0.5","3.6","0.9","0");
#08
#my @data=("0","190.9","116.7","62.6","47.4","33.2","13.9","12.9","18.3","14.4","9.7","7.9","3.7","6.2","3.0","2.8","8.3","1.6","5.7","2.5","0.5","0","1.0","0","2.0","0","0.7");
#09
#my @data=("0","202.7","119.6","66.4","39.3","36.1","13.7","10.4","20.7","13.6","10.3","7.0","5.1","9.0","1.1","6.3","5.4","3.9","4.5","0","3.6","4.1","1.1","0.6","1.3","1.9","1.3");
#10
#my @data=("0","195.7","124.5","70.1","51.5","44.7","14.1","18.2","21.3","16.0","8.1","9.7","6.4","10.1","2.3","5.3","4.6","1.9","5.4","0","1.4","2.3","2.2","0","2.1","1.5","0.6");
#11
#my @data=("0","185.2","122.0","66.0","48.0","36.9","13.3","21.3","20.3","13.0","18.3","9.1","6.3","8.4","7.0","1.7","2.9","0.7","4.5","0.9","1.8","4.9","1.7","2.0","1.9","0.7","0");
#12
#my @data=("0","179.8","96.5","61.9","50.8","36.7","9.1","20.1","18.3","9.8","10.0","10.3","6.6","8.5","3.6","2.1","7.0","3.8","6.6","3.7","0","0","3.1","0.7","3.1","1.7","0");
#13
#my @data=("0","194.4","116.1","59.8","41.6","29.0","20.6","12.8","11.9","13.2","14.1","13.0","1.2","6.1","3.1","5.4","3.6","6.1","5.3","0","1.9","0.7","0","0","2.5","1.6","1.5");
#14
#my @data=("0","191.4","95.8","63.0","50.7","36.9","17.2","17.2","17.3","3.6","0.8","8.0","4.2","2.3","5.4","2.2","2.4","1.8","3.9","1.5","0","2.3","0.4","1.0","0","1.0","0.2");
#15
#my @data=("0","178.7","125.2","62.3","43.3","31.1","25.5","15.9","20.5","13.3","9.9","8.6","5.4","12.3","8.1","1.0","6.9","0.3","4.9","4.0","0.9","0.7","0.9","0.4","2.1","0.4","1.3");
#16
#my @data=("0","179.4","112.1","67.7","34.9","36.3","11.0","7.2","11.3","15.6","4.8","10.1","7.7","8.4","0","6.3","7.2","3.2","6.6","1.4","0","1.0","3.4","0","2.4","2.1","1.0");
#17
#my @data=("0","176.7","106.8","64.5","44.3","37.9","15.7","18.8","24.5","9.3","6.7","8.5","3.3","2.6","2.0","0.9","4.4","3.4","6.0","1.6","0","1.5","1.9","0.9","1.1","0","1.7");
#18
#my @data=("0","171.2","108.1","66.2","42.3","27.7","18.9","22.1","9.5","9.3","7.7","11.3","7.5","12.6","3.8","1.5","6.2","4.1","4.3","0","1.9","0","4.7","0.2","2.8","0","1.7");
#19
#my @data=("0","169.9","114.4","66.6","41.2","25.5","21.3","18.9","17.8","9.7","5.3","9.1","5.1","5.3","3.7","1.4","5.4","1.8","1.5","0","0.5","2.3","2.5","1.2","0.9","2.1","0");
#20
#my @data=("0","193.6","123.3","51.5","40.1","33.1","20.2","23.8","19.6","15.2","9.3","5.9","6.3","7.7","0.5","1.9","5.2","1.7","2.7","1.8","3.5","0.9","0","2.0","0","1.9","1.5");
#21
#my @data=("0","203.6","132.6","73.6","42.3","44.6","11.3","8.5","16.1","18.7","7.9","6.7","6.4","10.0","5.0","6.0","5.9","3.0","9.5","0","0","4.7","0","2.0","0.8","0.2","0.7");
#22
#my @data=("0","190.8","105.9","54.8","42.8","25.7","22.8","14.5","15.0","15.3","3.2","10.9","6.9","5.2","3.4","6.5","4.8","1.6","6.0","0.1","1.7","1.1","0.7","3.8","3.2","0.3","0.8");
#23
#my @data=("0","184.7","100.5","65.5","58.4","37.8","19.0","19.0","19.8","14.2","4.8","13.0","7.6","5.3","8.9","3.7","4.9","0.4","6.2","0.6","4.0","0.5","0.4","2.6","1.3","1.9","2.1");
#24
#my @data=("0","162.9","115.2","78.7","43.9","38.2","19.7","9.5","21.5","6.6","3.9","5.0","8.8","6.3","2.2","3.7","5.1","0.4","4.8","1.4","4.1","0","2.2","0.1","1.5","0.1","0");
#25
#my @data=("0","194.0","129.7","55.7","53.0","35.1","22.8","11.2","18.2","15.2","10.4","10.2","6.1","10.2","4.8","2.2","9.9","2.6","3.1","0","0.6","0.6","0.4","3.7","3.2","0","1.2");
#26
#my @data=("0","178.7","94.5","51.3","34.2","34.8","20.9","22.3","22.8","12.9","8.1","13.2","2.8","6.6","2.0","5.5","7.4","1.5","5.3","3.6","0.2","2.6","3.0","0.7","2.5","1.1","1.1");
#27
#my @data=("0","178.1","114.4","64.0","50.6","34.8","19.7","15.2","17.5","7.9","7.2","8.8","4.6","8.3","5.0","1.6","8.6","1.4","4.4","4.2","0","2.8","0.7","2.2","2.7","2.2","0");
#28
#my @data=("0","189.0","107.7","57.3","46.5","33.1","18.3","20.8","19.3","16.2","10.0","9.4","4.4","7.8","6.4","3.7","3.4","0","6.6","2.8","0.9","2.6","1.9","1.9","1.8","2.0","0.2");
#29
#my @data=("0","189.1","102.1","69.9","39.7","46.2","17.8","15.2","17.6","15.1","10.4","13.5","1.0","6.5","7.1","3.4","6.0","0.9","5.8","0.7","0","0.3","1.2","0.7","3.6","1.9","0.3");
#30
#my @data=("0","192.1","131.3","64.2","46.3","34.8","15.8","20.7","24.2","10.6","5.4","11.1","5.7","8.4","2.1","5.0","8.6","2.0","7.0","0","2.4","0.8","0","2.4","2.8","0.4","0.6");
#31
#my @data=("0","175.7","106.0","67.9","52.7","41.6","26.6","10.3","20.0","15.4","9.1","4.1","3.9","6.6","4.4","5.2","8.3","2.6","5.1","0.2","1.4","3.1","0.3","0.3","3.7","0","1.9");
#32
#my @data=("0","193.4","109.7","48.4","55.5","38.1","9.3","15.2","19.5","12.1","6.9","12.2","6.8","8.7","6.3","2.4","5.8","0","7.1","0","0.5","0","0.8","0.8","2.2","2.1","2.2");
#33
#my @data=("0","188.0","105.5","58.1","66.3","27.1","25.4","13.5","19.6","6.0","0","11.3","1.6","8.3","5.6","3.6","2.1","0","4.9","1.5","1.1","2.6","0","0.8","1.9","1.9","2.1");
#34
#my @data=("0","174.4","116.3","62.7","43.7","33.6","16.3","16.1","15.1","13.5","8.1","9.6","5.1","8.7","6.2","3.0","5.5","1.0","6.1","1.0","0","1.3","0.7","0","2.0","0.7","0.8");
#35
#my @data=("0","188.0","117.7","59.0","49.0","41.3","18.0","13.6","16.1","15.8","13.7","8.9","5.2","5.9","4.2","5.0","4.5","2.4","4.7","4.5","2.2","2.9","0.7","0","2.9","0.9","1.1");
#36
#my @data=("0","189.8","101.7","71.3","49.9","37.4","17.5","20.1","15.3","14.7","11.1","9.3","6.7","7.7","5.4","3.5","6.2","2.4","5.3","0","2.7","3.7","0","3.1","1.0","0.6","1.0");
#37
#my @data=("0","204.3","110.5","63.8","54.6","33.9","14.5","16.0","20.8","14.6","8.7","10.4","0.2","5.0","3.5","0","8.9","1.6","7.5","1.9","0.7","3.3","2.1","0.8","2.2","1.2","1.7");
#38
#my @data=("0","194.2","105.8","61.2","53.3","33.3","18.3","16.0","17.2","14.2","4.0","10.9","3.9","10.9","0.8","0","5.2","3.9","3.4","2.0","0.1","1.2","1.2","0","2.2","0","0.6");
#39
#my @data=("0","170.2","123.4","61.7","57.9","27.8","17.6","18.4","20.9","3.4","10.9","4.1","5.2","4.4","5.9","4.2","5.2","2.8","6.3","1.6","1.8","1.7","1.5","0","0.2","0.6","1.2");
#40
#my @data=("0","188.1","111.7","52.0","42.7","31.3","18.4","7.0","18.8","16.1","3.3","10.5","1.4","9.9","3.6","3.3","6.9","0.6","5.2","0","3.2","2.1","0.1","1.9","2.3","0.2","0.7");
#41
#my @data=("0","177.1","119.5","69.1","51.5","30.6","17.4","14.5","26.6","10.5","6.4","8.7","0.5","6.0","1.5","2.2","5.6","1.4","4.9","0","1.8","3.5","0.8","1.4","2.5","2.0","0.6");
#42
#my @data=("0","186.5","122.3","60.6","40.5","33.2","10.7","20.4","17.7","14.0","8.6","6.5","0.9","5.1","5.3","4.0","8.0","0.2","4.3","2.1","1.5","3.7","0","0","4.0","2.0","0.8");
#43
#my @data=("0","201.1","110.9","64.5","52.3","35.7","17.8","14.2","21.4","12.5","7.2","2.6","9.2","7.6","4.2","5.1","7.6","2.6","2.2","0","0.8","2.6","1.3","0.5","3.7","2.4","0.5");
#44
#my @data=("0","204.5","116.7","62.1","47.0","29.9","21.2","10.5","16.2","17.1","10.9","12.3","7.6","6.4","3.1","3.6","3.4","2.4","6.3","2.6","0.3","1.8","0.5","1.3","2.1","1.0","1.0");
#45
#my @data=("0","184.4","130.9","61.9","48.1","38.3","17.2","20.4","14.7","9.3","10.4","7.5","4.4","5.1","3.9","5.5","6.4","3.5","3.4","3.0","2.6","3.9","0","3.9","2.0","0.8","2.0");
#46
#my @data=("0","177.8","124.4","77.3","45.8","34.6","19.6","13.2","14.6","13.1","7.9","7.8","4.9","7.1","2.3","5.3","4.2","0","5.8","0.7","1.5","1.4","1.3","3.2","2.4","3.6","1.6");
#47
#my @data=("0","178.8","85.0","58.9","61.2","36.1","20.6","16.2","13.0","15.5","9.4","11.4","3.8","10.7","3.1","7.7","6.4","3.8","4.8","2.3","0","1.0","4.0","0","3.5","0.7","1.5");
#48
#my @data=("0","213.0","106.0","57.5","50.4","31.5","11.1","19.1","17.2","13.0","6.9","9.3","3.5","10.0","6.8","3.4","5.7","2.3","4.7","0.9","1.4","1.5","0","0.7","2.0","1.9","1.3");
#49
#my @data=("0","198.0","119.6","50.2","48.7","31.1","21.7","9.9","15.2","13.7","7.0","9.6","6.8","6.9","0","3.4","4.9","1.4","6.2","0","0","3.3","0.1","1.7","4.1","1.6","3.2");
#50
#my @data=("0","172.6","118.0","72.1","47.9","36.2","23.8","8.9","19.9","8.9","9.2","7.3","3.4","7.3","3.4","10.6","7.2","1.4","4.3","2.7","1.9","2.2","0","0.1","2.3","1.5","0.6");
#51
#my @data=("0","203.0","102.4","64.2","51.6","30.3","21.4","16.4","17.8","12.9","6.6","9.4","8.1","10.2","4.9","6.2","2.4","5.5","2.9","0.5","1.8","2.0","2.0","1.4","2.8","0.3","1.7");
#52
#my @data=("0","165.8","119.3","66.2","42.1","30.3","17.4","14.7","24.9","11.3","9.7","7.5","0","10.2","5.5","4.7","6.7","1.6","2.9","0","0","4.0","1.5","1.7","3.4","1.4","1.5");
#53
#my @data=("0","188.9","104.6","65.6","43.3","27.8","12.8","15.2","15.4","8.0","7.6","13.1","6.5","4.7","9.6","2.2","7.4","3.1","5.0","0.9","0","2.9","0.9","1.6","2.4","2.4","1.9");
#54
#my @data=("0","177.0","123.5","75.5","55.4","32.3","19.4","17.2","10.7","11.4","9.6","9.6","1.1","6.1","0","4.5","6.1","2.6","2.8","1.6","1.1","2.1","0.1","0","2.1","0.2","1.3");
#55
#my @data=("0","180.7","120.0","68.4","54.4","32.5","14.7","19.0","12.8","8.3","7.8","13.4","0","6.5","1.9","0","6.3","0.5","6.3","0","1.6","1.2","3.6","1.7","3.4","0","2.7");
#56
#my @data=("0","188.2","108.4","64.4","40.7","38.8","21.3","12.0","20.9","16.5","6.8","9.9","5.7","7.1","3.6","7.7","5.0","0.6","7.3","2.1","0","3.8","1.3","0.5","2.2","2.1","1.5");
#57
#my @data=("0","178.6","117.5","58.9","51.2","37.7","26.9","11.9","22.1","13.2","17.7","10.8","5.0","8.5","6.0","6.4","7.5","2.1","6.6","2.9","0","3.0","0.9","1.4","0.6","1.8","0");
#58
#my @data=("0","178.7","116.0","61.2","49.1","25.2","18.1","13.2","22.2","14.4","9.5","11.7","2.4","4.8","1.2","4.1","1.0","3.6","5.9","0.4","0","0","0.6","1.1","2.2","1.2","1.7");
#59
#my @data=("0","165.4","102.9","71.1","45.4","33.9","20.5","14.9","19.0","5.0","5.2","9.9","4.4","7.9","2.5","3.1","6.7","4.8","6.5","1.6","1.3","3.3","1.9","0.2","0.7","0.8","0");
#60
#my @data=("0","176.4","109.3","62.1","51.6","27.3","16.2","16.3","15.7","10.0","9.5","11.7","5.7","6.3","3.9","3.0","8.0","2.4","6.6","0.2","0.2","1.2","2.6","1.9","2.7","1.3","0.4");
#61
#my @data=("0","187.7","115.4","56.1","54.0","29.6","15.1","10.9","20.2","21.8","8.0","12.6","3.9","9.0","4.0","4.8","6.1","0.7","5.4","0.3","1.7","2.9","0","0.4","2.4","1.8","1.2");
#62
#my @data=("0","210.6","86.0","51.0","50.9","30.9","12.6","19.2","15.3","6.7","11.7","10.8","6.5","7.5","2.7","3.2","6.8","2.9","5.1","1.3","1.1","2.4","2.1","0","0","0.5","0.9");
#63
#my @data=("0","181.7","120.4","57.4","50.2","27.0","17.9","20.2","20.5","13.6","3.7","10.2","3.3","4.3","2.0","3.5","8.0","4.0","5.6","0","2.2","0.2","1.2","0","1.5","2.2","2.9");
#64
#my @data=("0","177.8","96.0","55.7","42.3","37.8","26.7","14.0","21.4","13.8","10.3","9.6","3.9","9.7","5.1","6.7","4.7","1.9","5.7","1.8","0.6","1.6","0.3","0","2.7","0","0");
#65
#my @data=("0","186.2","97.6","62.5","60.3","31.4","17.3","15.8","22.6","11.3","14.1","13.1","1.7","4.3","4.6","3.6","6.9","4.6","4.1","1.1","0","2.5","0.5","0","4.6","1.1","2.3");
#66
#my @data=("0","178.8","104.2","52.3","60.7","34.0","19.0","23.9","23.1","12.7","4.6","8.2","2.5","4.0","5.0","4.4","6.1","2.2","3.5","5.0","0","2.2","2.9","0.7","0.7","1.1","1.5");
#67
#my @data=("0","167.9","115.3","58.2","48.9","35.2","22.9","25.6","16.5","12.7","3.1","6.8","5.7","8.6","1.5","8.9","6.6","2.9","4.7","0","1.3","0","1.5","2.1","2.8","0","1.3");
#68
#my @data=("0","176.3","117.4","44.4","45.8","43.5","18.2","23.1","22.7","17.4","7.3","7.8","2.6","7.1","4.7","3.3","3.4","1.8","2.4","3.1","3.7","0","0","0.1","1.6","1.8","0.8");
#69
#my @data=("0","179.0","96.1","61.9","58.5","31.1","14.3","12.4","27.6","10.3","6.7","10.5","0.9","4.7","3.2","3.8","2.7","2.5","6.0","0","0","2.9","1.2","0.5","2.2","0","0.4");
#70
#my @data=("0","188.0","113.4","59.2","45.2","30.7","11.9","8.6","25.8","10.3","8.9","5.7","0.8","7.3","5.4","5.2","7.0","0.3","7.4","2.0","0.9","2.5","3.1","0.7","2.3","0.6","0.3");
#71
#my @data=("0","191.4","113.9","46.4","54.7","35.7","20.0","19.7","18.6","11.5","4.5","9.5","8.5","7.5","7.6","7.3","2.3","1.6","4.8","0","1.6","1.9","1.1","1.9","1.5","0.7","1.1");
#72
#my @data=("0","186.7","117.6","65.8","60.9","42.0","10.4","27.0","15.3","9.2","5.8","14.5","4.5","6.6","5.0","0.7","6.1","3.6","4.2","0.5","2.4","3.3","0","3.1","1.1","0.5","0.7");
#73
#my @data=("0","183.8","101.8","61.7","53.0","40.4","20.9","18.4","23.7","8.3","6.9","8.1","5.5","5.1","1.3","0.8","6.9","5.1","5.3","1.8","1.1","2.5","1.2","0.4","1.5","0","0");
#74
#my @data=("0","151.9","114.0","58.9","55.4","34.8","23.4","15.2","19.0","9.7","10.5","11.3","6.6","6.5","0.4","2.7","6.9","1.1","5.3","3.0","0","4.6","0","2.5","1.0","1.7","1.3");
#75
#my @data=("0","200.4","112.8","72.5","47.8","29.9","12.2","17.2","22.3","15.0","10.5","10.8","4.9","5.3","3.4","2.5","7.4","4.2","5.9","0.6","2.8","2.1","1.2","1.3","3.6","0.8","1.8");
#76
#my @data=("0","185.2","143.7","57.6","56.8","29.3","11.7","15.3","20.9","10.0","6.8","11.9","6.6","9.6","4.6","3.3","5.1","1.8","3.2","0.6","0.8","0.5","2.2","1.0","1.2","1.4","2.3");
#77
#my @data=("0","177.3","120.0","63.6","39.3","31.2","17.4","18.0","17.8","5.9","13.2","10.7","3.2","5.3","4.7","2.9","8.6","0.4","2.8","3.3","2.4","3.2","0.5","0","2.5","1.4","0.2");
#78
#my @data=("0","186.2","101.6","57.0","50.4","29.1","30.2","15.5","17.7","13.1","7.3","7.5","2.3","6.0","0.8","4.6","8.8","2.7","2.4","2.4","0.5","2.4","0.5","1.2","2.3","1.5","1.4");
#79
#my @data=("0","192.9","123.0","76.6","54.2","30.4","15.8","11.3","19.8","13.0","10.0","10.7","0.6","4.1","1.8","2.7","7.2","0","4.5","2.1","4.7","0.5","2.1","4.3","5.1","1.1","0.7");
#80
#my @data=("0","207.5","106.1","59.1","52.5","29.9","15.4","23.2","20.5","14.5","1.2","8.2","0","5.4","2.7","7.5","1.2","3.5","7.1","0","2.0","3.4","0","1.8","2.9","1.4","0.2");
#81
#my @data=("0","175.6","99.8","53.1","42.7","36.3","20.5","13.2","19.4","13.9","13.0","9.0","3.4","2.2","2.2","5.3","5.5","3.3","3.0","0.8","2.3","2.4","2.2","1.4","4.1","0","0.1");
#82
#my @data=("0","172.0","112.3","53.7","37.3","31.0","11.1","11.5","16.7","13.1","6.3","11.6","7.6","3.6","2.3","3.6","5.4","3.7","7.2","2.0","0","2.4","0.3","1.2","1.5","2.8","0.8");
#83
#my @data=("0","206.6","99.9","54.5","57.6","40.3","12.2","19.3","18.2","9.2","7.8","10.5","6.1","7.8","2.9","6.0","5.2","0","5.2","0.9","0.2","3.6","3.1","0","2.1","1.3","0");
#84
#my @data=("0","170.8","133.8","62.6","53.2","37.7","18.8","14.6","18.8","17.6","15.7","6.4","5.7","5.9","3.1","3.6","5.6","2.5","5.9","0","0","4.5","0.5","0.5","2.2","1.3","0.8");
#85
#my @data=("0","191.8","118.7","64.3","47.8","35.3","13.0","25.6","21.0","10.4","5.3","18.1","4.4","5.8","6.4","0","6.4","2.0","6.9","2.8","1.3","0.8","1.0","3.8","1.0","1.5","2.2");
#86
#my @data=("0","185.7","109.3","56.0","53.5","43.2","24.6","16.3","14.4","11.3","7.1","14.0","8.7","7.6","4.0","3.2","4.8","5.1","4.1","0","2.3","0.5","3.1","2.6","2.0","1.7","2.2");
#87
#my @data=("0","195.5","101.4","71.5","54.4","31.9","22.8","15.8","16.3","8.2","7.0","12.8","5.9","10.1","5.0","7.5","2.4","3.5","3.4","2.4","0","3.3","0","1.0","1.3","0.2","2.0");
#88
#my @data=("0","181.5","92.6","74.2","55.9","40.3","14.7","15.0","29.6","16.2","5.2","9.0","2.4","7.5","2.0","4.1","3.8","1.9","3.7","0.8","1.8","2.0","1.4","0.3","0.1","1.2","3.0");
#89
#my @data=("0","166.6","106.8","68.5","55.6","36.7","16.1","22.7","15.3","17.2","8.5","3.0","9.0","8.4","4.9","3.3","4.8","1.9","2.7","0","1.1","2.1","0.5","0.3","1.5","0","0");
#90
#my @data=("0","189.1","103.3","67.4","49.7","33.5","17.4","14.4","19.4","14.7","0","7.9","9.8","7.3","2.4","1.6","7.4","1.9","8.8","2.6","0.9","3.5","0","0","1.8","2.6","0");
#91
#my @data=("0","182.2","101.3","69.8","52.1","21.8","20.2","15.1","25.7","7.4","3.8","10.2","3.7","11.7","6.0","1.2","6.9","3.6","3.7","0.1","0.7","1.7","0","2.3","2.4","0.9","1.4");
#92
#my @data=("0","192.2","115.9","62.6","45.4","37.8","21.5","23.4","20.6","9.2","7.0","7.4","8.3","6.8","2.6","3.4","3.1","4.5","4.7","0","0.6","3.9","1.2","1.2","1.6","0","1.1");
#93
#my @data=("0","190.8","119.0","61.4","61.4","42.7","21.1","14.6","20.3","16.6","6.0","10.1","0","9.8","4.7","3.0","0.7","2.7","2.3","2.8","0.1","3.2","0.7","0.8","2.5","2.0","0");
#94
#my @data=("0","175.7","106.5","69.8","46.9","29.8","15.6","9.3","20.3","11.3","14.4","12.4","1.1","5.3","1.4","1.3","7.1","3.1","6.1","2.4","1.1","2.9","0","1.1","2.7","1.6","1.6");
#95
#my @data=("0","193.7","107.8","66.4","49.2","31.1","19.2","11.2","21.1","14.8","13.8","9.6","5.1","6.4","5.4","1.5","7.5","0.8","6.0","1.6","1.0","2.3","0","0.9","2.4","1.5","0.8");
#96
#my @data=("0","182.3","105.2","56.9","48.9","39.1","14.0","13.3","18.4","12.7","5.8","7.9","4.0","10.7","6.2","4.0","6.8","7.0","5.6","0","1.2","2.6","0.3","2.5","1.1","0.5","2.8");
#97
#my @data=("0","200.2","97.2","52.5","48.0","41.8","20.8","17.8","20.4","8.3","5.8","5.6","0.6","8.1","6.6","6.6","5.9","4.8","3.8","0","1.8","3.4","0","0.6","3.0","1.8","0.8");
#98
#my @data=("0","188.6","114.4","63.5","45.1","33.4","13.8","21.4","17.8","14.6","5.6","10.7","10.4","4.5","0.9","4.7","2.9","0.8","3.6","0","0","2.1","2.2","0","1.5","0.2","1.3");
#99
#my @data=("0","181.9","120.6","66.1","45.2","30.0","20.8","22.8","12.0","10.6","9.4","8.0","3.8","8.5","8.2","0","4.9","0","7.9","0.7","2.1","2.2","0.2","0.6","2.8","1.7","2.0");


#2, synthetic contest set
#626 contacts, 111 calls, synthetic
#@data originally started at "19", the number of people who were worked 1 time
#the simulation determines the number of people who were worked 0 times (which is 0 here)
#my $contacts = 626;
#my @data=("0","19", "16", "13", "9", "9", "4", "5", "7", "0", "2", "3", "2", "1", "1", "1", "0", "1", "3", "1", "1", "1", "1", "0", "0", "1", "0", "0", "1");

#monte carlo data
#my @data=("0.2", "19.6", "16.5", "12.9", "9.1", "9.1", "4.2", "4.8", "6.8", "-0.0", "1.9", "3.1", "2.1", "0.9", "1.1", "1.0", "0.1", "1.1", "3.1", "0.9", "1.1", "0.8", "0.9", "-0.0", "-0.0", "1.0");

#my @data=("-0.0", "18.5", "16.6", "13.1", "8.3", "8.5", "4.4", "5.0", "7.1", "-0.0", "2.1", "2.9", "2.0", "1.1", "0.9", "1.0", "-0.0", "0.9", "3.2", "1.0", "0.9", "1.1", "1.1", "0.2", "0.2", "1.0");

#my @data=("0.7", "19.2", "16.3", "13.0", "8.5", "8.5", "4.1", "4.9", "6.9", "0.0", "2.0", "2.9", "1.9", "0.9", "1.1", "0.8", "0.0", "1.0", "2.9", "0.9", "1.0", "1.0", "1.0", "-0.0", "-0.0", "1.0");

#my @data=("-0.0", "19.2", "15.8", "13.1", "9.0", "8.8", "4.1", "4.8", "6.8", "0.0", "2.0", "3.0", "1.9", "1.0", "1.0", "1.0", "-0.0", "1.2", "3.0", "1.1", "1.2", "1.0", "1.2", "0.0", "-0.0", "1.0");

#my @data=(" 0.3", "18.7", "16.1", "13.2", "8.9", "8.7", "4.1", "4.9", "7.0", "0.2", "2.0", "2.9", "2.1", "1.0", "1.1", "1.0", "-0.0", "0.9", "2.9", "0.9", "0.9", "1.0", "1.0", "0.1", "0.1", "1.0");

#my @data=("-0.0","18.8","15.8","12.7","9.2","8.5","3.9","4.9","7.4","0.1","1.9","3.0","1.9","1.0","0.9","1.1","-0.0","1.0","3.1","1.1","1.1","1.0","0.9","0.0","-0.0","1.0");

#for ($i=0;$i < $#data; $i++)
#	{
#	print "@data[$i] \n";
#	}

#depending on the model, uncomment one of these 
my $population = @ARGV[1];	#probability is a function
#my $population = @ARGV[10];	#probability is determined from 10 bins

my $member_number;
my $i;
my $end_loop;
my ($accum_error,$accum_error_15);
my $expn = 2.7182818;

my $QSOs_made = 0;
my $QSOs_missed = 0;
my $accum_QSOs_missed = 0;
my $operators_seen=0;
my $lowest_common_max_contacts=$contacts; #at the end of the simulation will be the number seen in all repeats 
my $operators_seen_sum=0;
my $operators_seen_sum_of_squares=0;
my $loops_done = 0;
my $exponent;
my ($error, $error_numerator, $error_denominator);
my (@array_QSOs_made_by_member, @array_QSOs_missed_by_member); 
my (@array_QSOs_made_by_density, @array_QSOs_missed_by_density);

my $simulations=1024;
my $simulation;
my (@array_QSOs_made_by_density_sum, @array_QSOs_made_by_density_sum_of_squares); 
my (@array_QSOs_missed_by_density_sum, @array_QSOs_missed_by_density_sum_of_squares); 

sub initialise_arrays
	{
	#perl arrays start at 0. Data starts at 1
	#arrays in perl are initialised to null. 
	#if a particular operator is not worked, I want this to be initiallised to 0
	for ($i=0; $i <= $population; $i++)
		{
		@array_QSOs_made_by_member[$i]=0;
		@array_QSOs_made_by_density[$i]=0;
		@array_QSOs_missed_by_member[$i]=0;
		@array_QSOs_missed_by_density[$i]=0;
		#print "initialise_arrays: @array_QSOs_made_by_member[$i], @array_QSOs_missed_by_member[$i];";
		}
	}

sub offer_contacts
	{
	#for ($i=0; $i< $contacts; $i++)
	$QSOs_made=0;
	$QSOs_missed=0;
	$loops_done=0;
	#printf ("offer_contacts start: QSOs_made=%6d contacts= %6d \n", $QSOs_made, $contacts);
	while ($QSOs_made < $contacts)
		{
		#pick one of the population
		$member_number = int (rand ($population)); #goes from 0..(popln-1)
	
		#do we count it?
		# case 2: chance of QSO being made is proportional to an exp function
		$exponent = @ARGV[1]*@ARGV[0]*$member_number/$population;
		$exponent = sqrt($exponent);
		$exponent = -$exponent;
		#printf ("offer_contacts: member %4d ", $member_number);
		if ($expn**($exponent) > rand(1))
			{
			@array_QSOs_made_by_member[$member_number]++;
			$QSOs_made++;
			#printf ("made \n");
			}
		else	
			{
			@array_QSOs_missed_by_member[$member_number]++;
			$QSOs_missed++;
			$accum_QSOs_missed++;
			#printf ("missed \n");
			}
		$loops_done++;
		}
	#printf ("offer_contacts exit:QSOs_made=%6d QSOs_missed= %6d loops_done=%6d\n", $QSOs_made, $QSOs_missed, $loops_done);
	}

sub calculate_frequency_distribution
	{
	#how many members were seen/made QSOs 1,2,3..n times?
	$operators_seen=0;	# if we find a ham that wasn't worked (@array_QSOs_made_by_member[$i] == 0),
					# then we decrement $operators_seen.
	for ($i=0;$i < $population;$i++)
		{
		#print "calculate_frequency_distribution: $i, @array_QSOs_made_by_member[$i], @array_QSOs_missed_by_member[$i]";
		$QSOs_made += @array_QSOs_made_by_member[$i];
		$QSOs_missed += @array_QSOs_missed_by_member[$i];
		#if this member $i==10 was seen 3 times, then increment @array_QSOs_made_by_density[3]
		#all we know now is how many people were worked 3 times.
		#we loose the identity of the members who were worked 3 times.
		@array_QSOs_made_by_density[@array_QSOs_made_by_member[$i]]++;
	
		#this is a bit trickier
		#if member 10 missed 100 qsos, but made 3 QSOs, 
		#then we want to record that the people who were worked 3 times 
		#missed 100 qsos
		@array_QSOs_missed_by_density[@array_QSOs_made_by_member[$i]] += @array_QSOs_missed_by_member[$i];
	
		if (@array_QSOs_made_by_member[$i] != 0) 
			{
			$operators_seen++;
			#print "$operators_seen ";
			}
		#print "\n";
		}
	#print "calculate_frequency_distribution:QSOs_made= $QSOs_made QSOs_missed= $QSOs_missed operators_seen= $operators_seen \n";
	$operators_seen_sum += $operators_seen;
	$operators_seen_sum_of_squares += $operators_seen*$operators_seen;
	}

sub print_output
	{
	#output
	$accum_error=0;
	$operators_seen=0;
	#print to the largest valid data in either @data or @array_QSOs_made_by_density
	# $#data is the maximum times any operator was worked in the contest
	#the @array_QSOs_made_by_density has been initialised to 0 so can't
	#use $#array_QSOs_made_by_density
	#count backwards till find non-zero element
	
	$end_loop=$population;
	while (@array_QSOs_made_by_density[$end_loop] == 0)
		{
		$end_loop--;
		}
	$end_loop++;
	#set the high water mark for freuqency distribution
	if ($end_loop < $lowest_common_max_contacts) {$lowest_common_max_contacts = $end_loop;}
	if ($#data > $end_loop) {$end_loop = $#data;}
	for ($i=0;$i < $end_loop;$i++)
		{
		#calculate error fn
		#printf ("%3d ",@data[$i]);
		#printf ("%3d ",@array_QSOs_made_by_density[$i]);
		if (@data[$i] == 0) { $error_denominator = 1; }
		else { $error_denominator = @data[$i];}
		if (@array_QSOs_made_by_density[$i] == 0){$error_numerator = 1;}
		else { $error_numerator = @array_QSOs_made_by_density[$i];}
		$error = $error_numerator/$error_denominator;
		if ($error < 1) {$error = 1/$error;}
		if ((@data[$i] == 0) &&  (@array_QSOs_made_by_density[$i] == 0))
			{;}
		else 	{
			#ignore error for people being worked 0 times
			if ($i != 0)
				{
				#$accum_error += $error*$error;
				$accum_error += $error*@data[$i];
				}
			}
		#printf ("output %3d %4d %4d %4d %6.1f ",$i, @data[$i], @array_QSOs_made_by_density[$i], @array_QSOs_missed_by_density[$i], $error);
	
		#don't count the operators who were seen zero times
		if ($i > 0) {$operators_seen += @array_QSOs_made_by_density[$i];}
		#printf ("%4d ", $operators_seen);
		#printf ("%6.1f ", $accum_error) ;

		#accumulate for multiple simulations
		
		@array_QSOs_made_by_density_sum[$i] += @array_QSOs_made_by_density[$i];
		@array_QSOs_made_by_density_sum_of_squares[$i] += @array_QSOs_made_by_density[$i] 
								* @array_QSOs_made_by_density[$i];
		@array_QSOs_missed_by_density_sum[$i] += @array_QSOs_missed_by_density[$i];
		@array_QSOs_missed_by_density_sum_of_squares[$i] += @array_QSOs_missed_by_density[$i] 
								* @array_QSOs_missed_by_density[$i];
		#printf ("sums %4d %8d %4d %8d", @array_QSOs_made_by_density_sum[$i], 
		#				@array_QSOs_made_by_density_sum_of_squares[$i], 
		#				@array_QSOs_missed_by_density_sum[$i], 
		#				@array_QSOs_missed_by_density_sum_of_squares[$i]);
		#print "\n";
	
		$loops_done++;
		}
	#print "QSO's_made $QSOs_made loops_done $loops_done operators_seen $operators_seen \n";
	}

sub calculate_variances
	{
	my ($QSOs_made_mean, $QSOs_missed_mean, $operators_seen_mean);
	my ($QSOs_made_variance, $QSOs_missed_variance, $operators_seen_variance);
	my ($QSOs_made_standard_deviation, $QSOs_missed_standard_deviation, $operators_seen_standard_deviation);
	my ($QSOs_made_average_square_of_sums, $QSOs_missed_average_square_of_sums, $operators_seen_average_square_of_sums);
	my ($QSOs_made_average_sum_of_squares, $QSOs_missed_average_sum_of_squares, $operators_seen_average_sum_of_squares);
	my ($miss_to_made_ratio, $miss_to_made_ratio_standard_deviation);

	$accum_error = 0;

	printf ("calculate_variances:contacts data calc_made +/- calc_missed +/- error acc_er miss/made +/- \n");
	for ($i=0;$i < $lowest_common_max_contacts;$i++)
		{
		#accumulate for multiple simulations
		#calculate mean
		$QSOs_made_mean = @array_QSOs_made_by_density_sum[$i]/$simulations; 
		#$QSOs_made_average_sum_of_squares = @array_QSOs_made_by_density_sum_of_squares[$i]/($simulations-1); 
		$QSOs_made_average_sum_of_squares = @array_QSOs_made_by_density_sum_of_squares[$i]/$simulations; 
		$QSOs_made_average_square_of_sums = @array_QSOs_made_by_density_sum[$i] 
						* @array_QSOs_made_by_density_sum[$i]/($simulations*$simulations); 
		$QSOs_made_variance = $QSOs_made_average_sum_of_squares - $QSOs_made_average_square_of_sums;
		$QSOs_made_standard_deviation = sqrt($QSOs_made_variance);

		$QSOs_missed_mean = @array_QSOs_missed_by_density_sum[$i]/$simulations; 
		#$QSOs_missed_average_sum_of_squares = @array_QSOs_missed_by_density_sum_of_squares[$i]/($simulations-1); 
		$QSOs_missed_average_sum_of_squares = @array_QSOs_missed_by_density_sum_of_squares[$i]/$simulations; 
		$QSOs_missed_average_square_of_sums = @array_QSOs_missed_by_density_sum[$i] 
 						* @array_QSOs_missed_by_density_sum[$i]/($simulations*$simulations); 
		$QSOs_missed_variance = $QSOs_missed_average_sum_of_squares - $QSOs_missed_average_square_of_sums;
		$QSOs_missed_standard_deviation = sqrt($QSOs_missed_variance);
	
		#calculate error fn
		#printf ("%3d ",@data[$i]);
		#printf ("%3d ",$QSOs_made_mean);
		if (@data[$i] == 0) { $error_denominator = 1; }
		else { $error_denominator = @data[$i];}
		if ($QSOs_made_mean == 0){$error_numerator = 1;}
		else { $error_numerator = $QSOs_made_mean;}
		$error = $error_numerator/$error_denominator;
		if ($error < 1) {$error = 1/$error;}
		if ((@data[$i] == 0) &&  ($QSOs_made_mean == 0))
			{;}
		else 	{
			#ignore error for people being worked 0 times
			if ($i != 0)
				{
				# $accum_error += $error*$error;
				$accum_error += $error*@data[$i];
				}
			if ($i == 15)
				{
				#depending on input, output length and 
				#hence accum_error is different.
				#compare all at same place
				$accum_error_15 = $accum_error;
				}
			}
		#printf ("%8d %6.1f %6d %6.1f %6d", $i, $QSOs_made_mean, $QSOs_made_standard_deviation, $QSOs_missed_mean, $QSOs_missed_standard_deviation);
		printf ("calculate_variances: %6d %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f ", $i, @data[$i], $QSOs_made_mean, $QSOs_made_standard_deviation, $QSOs_missed_mean, $QSOs_missed_standard_deviation, $error, $accum_error );
		
		if ($i == 0) {print "     0      0\n";}
		else 	{
			$miss_to_made_ratio = $QSOs_missed_mean/($i * $QSOs_made_mean);
			$miss_to_made_ratio_standard_deviation = $miss_to_made_ratio*sqrt($QSOs_missed_variance/($QSOs_missed_mean*$QSOs_missed_mean) + $QSOs_made_variance/($QSOs_made_mean*$QSOs_made_mean));

			printf ("%6.3f %6.3f \n", $miss_to_made_ratio, $miss_to_made_ratio_standard_deviation);
			}
			
		}
	print "calculate_variances: input params @ARGV[0] @ARGV[1] \n";
	$operators_seen_mean = $operators_seen_sum/$simulations;
	$operators_seen_average_square_of_sums = $operators_seen_sum*$operators_seen_sum/($simulations*$simulations);
	#$operators_seen_average_sum_of_squares = $operators_seen_sum_of_squares/($simulations-1);
	$operators_seen_average_sum_of_squares = $operators_seen_sum_of_squares/$simulations;
	$operators_seen_variance = $operators_seen_average_sum_of_squares - $operators_seen_average_square_of_sums;	
	#print "$operators_seen_mean, $operators_seen_sum_of_squares, $operators_seen_square_of_sums";	
	printf ("calculate_variances: operators_seen = %6.1f +/- %6.1f QSOs missed = %6.1f \n", $operators_seen_mean, $operators_seen_standard_deviation, $accum_QSOs_missed/$simulations );
	printf ("calculate_variances: accum_error_15 %6.1f \n", $accum_error_15);

	printf ("output: %9.6f %6d seen=%6.1f +/- %3.1f missed=%6.1f acccum_error_15 %6.1f\n", @ARGV[0], @ARGV[1], $operators_seen_mean, $operators_seen_standard_deviation, $accum_QSOs_missed/$simulations, $accum_error_15 );
	print "\n";
	}

#-----------------------------------------------------------------------
#main:

#$#ARGV is the (number of args - 1), eg $foo arg1 arg2 gives $#ARGV==1
#(programming perl p138)

if (!($#ARGV == 1) )
	{
	print "simulate.pl (C) Joseph Mack NA3T 2000, jmack\@wm7d.net, released under GPL\n";
	print "usage:\n";
	print "simulate.pl scaling_factor population [-m]\n";
	print "scaling_factor - start with 0.02-0.03\n";
	print "population - expected number of operators\n";
	print "\n"; 
	print "Code has array for experimental data. Number of simulations is set in code too.\n"; 
	print "Calculates fit for input parameters\n";
	print "@ARGV[0], @ARGV[1], @ARGV[2], $#ARGV \n";
	exit;
	}


for ($simulation=0;$simulation<$simulations;$simulation++)
	{
	#print "\n";
	#printf "simulation $simulation \r";
	initialise_arrays;
	offer_contacts;
	calculate_frequency_distribution;
	print_output;
	}
calculate_variances;
#------------------------------------