RANDOMIZE(time)
close:cls: StartDate$=Today$: StartTime=time: StartTime$=time$
print "Confirm a directory (C:\Simulation results) for saving simulation results exists in your PC. The simulation produces .csv file for MS-EXCEL.":print
INPUT "Number of SubPopulations";NSP
input "Initial Number of Monandrous Female in each SubPopulation";INITNS
input "Initial Number of Polyandrous Female in each SubPopulation";INITNM
if InitNS>1 and InitNM>1 then MaxTrials=1000 else MaxTrials=50000
INPUT "Carrying Capacity of Each SubPopulation";KD
input "Mating Frequency of Polyandrous Female";MF
input "Maximum Fecundity of Females (>= 100)"; FDforMono
FDFORPOLYPERMATING=FDFORMONO/MF
input "Maximum Number of Calculating Generation";GN
input "population average of male quality (from 0 to 1)";AMQ
input "truncation point (from 0 to 0.5)";TRC
*MR
input "Migration Rate (from 0 to 0.9)";MR: if MR=1 then beep: goto *MR
input "% Costs of Polyandry (from 0 to 100)";cost
pMR=MR*100
F$=str$(MF)+"mating poly, male qual average"+str$(AMQ)+" trunc point"+str$(TRC)+", Poly"+str$(INITNM)+" vs Mono"+str$(INITNS)+","+str$(PMR)+"% dispers, cost"+str$(COST)+"%,"+str$(NSP)+"subpops,"+str$(MAXTRIALS)
open "C:\Simulation results\"+"continuousMQ, maxFD"+str$(FDFORMONO)+", "+F$+"iter.csv" for create as #1
print #1,"Number of SubPopulations= ";NSP;" / Mating Frequency of Polyandrous Female= ";MF;" / Costs of Polyandry= ";COST;" %";
print #1," / Migration Rate= ";MR;" / Potential for Population Increase (no cost cases)= ";FDFORMONO*AMQ;" times"
print #1, "Trial,"; "Initial Number of Monandrous Female in Each SubPopulation,"; "Initial Number of Polyandrous Female in Each SubPopulation,";
print #1, "Carrying Capacity of Each SubPopulation,"; "Mating Frequency of Polyandrous Female,"; "Number of Calculating Generation,"; "The Genaration at which Simulation ended,";
print #1, "Cumulative Number of Extinction of SubPopulation,";
print #1, "Cumulative Number of Monandry Extinction in SubPopulation,"; "Cumulative Number of Polyandry Extinction in SubPopulation,";
print #1, "Cumulative Number of Monandry Extinction in all SubPopulation,"; "Cumulative Number of Polyandry Extinction in all SubPopulation,";
print #1, "Cumulative Number of Both Extinct in all SubPopulation,"; "Cumulative Number of the Cases Polymorphism was kept at the End,";
print #1, "Total Number of Polyandrous Females in Metapopulation at the End,"; "Total Number of Monandrous Females in Metapopulation at the End"
DIM NS(GN+1,NSP+1),NM(GN+1,NSP+1),DispS(NSP+1),DispM(NSP+1),DISPS1(NSP+1),DISPM1(NSP+1),DISPS2(NSP+1),DISPM2(NSP+1)
TRIAL=1:SWIN=0:MWIN=0:BOTHEXTINCT=0:SSubpopExtinct=0:MSubpopExtinct=0:CumSSubpopExtinct=0:CumMSubpopExtinct=0:SubpopExtinct=0:CumSubpopExtinct=0
*STARTTRIALS
if INITNS=1 and INITNM>1 then for SUBPOP=1 to NSP: NS(1,SUBPOP)=0:NM(1,SUBPOP)=INITNM+1: next SUBPOP:NS(1,1)=1:NM(1,1)=initNM: T=1: SSUBPOPEXTINCT=0:MSUBPOPEXTINCT=0:SubpopExtinct=0: NMSUMATEND=0: NSSUMATEND=0: goto *STARTGEN
if INITNS>1 and INITNM=1 then for SUBPOP=1 to NSP: NS(1,SUBPOP)=INITNS+1:NM(1,SUBPOP)=0: next SUBPOP:NM(1,1)=1:NS(1,1)=INITNS: T=1: SSUBPOPEXTINCT=0:MSUBPOPEXTINCT=0:SubpopExtinct=0: NMSUMATEND=0: NSSUMATEND=0: goto *STARTGEN
for SUBPOP=1 to NSP: NS(1,SUBPOP)=INITNS:NM(1,SUBPOP)=INITNM: next SUBPOP: T=1: SSUBPOPEXTINCT=0:MSUBPOPEXTINCT=0:SubpopExtinct=0: NMSUMatEND=0: NSSUMatEND=0
*Startgen
for T=1 to GN
for Subpop=1 to NSP
SNUM=0:MNUM=0
if NS(T,SUBPOP)<1 then goto *nextNS
for I=1 to NS(T,SUBPOP)
*MQforMono
MQ=sqr(-2*log(rnd(1)))*cos(2*3.141592*rnd(1))/6+AMQ
if MQ1 then MQ=1
SNUM=SNUM+FDFORMONO*MQ
NEXT I
*nextNS
NS(T+1,SUBPOP)=SNUM
if NM(T,SUBPOP)<1 then goto *nextNM
for J=1 to NM(T,SUBPOP)
Mfitness=0
for MM=1 to MF
*MQforPoly
MQ=sqr(-2*log(rnd(1)))*cos(2*3.141592*rnd(1))/6+AMQ
if MQ1 then MQ=1
MFITNESS=MFITNESS+FDFORPOLYPERMating*MQ
next MM
MNUM=MNUM+MFITNESS
next J
*NEXTNM
NM(T+1,SUBPOP)=MNUM*(1-COST*0.01)
DISPS(SUBPOP)=int(MR*NS(T+1,SUBPOP)):NS(T+1,SUBPOP)=NS(T+1,SUBPOP)-DISPS(SUBPOP):DISPS1(SUBPOP)=int(DISPS(SUBPOP)/2):DISPS2(SUBPOP)=DISPS(SUBPOP)-DISPS1(SUBPOP)
DISPM(SUBPOP)=int(MR*NM(T+1,SUBPOP)):NM(T+1,SUBPOP)=NM(T+1,SUBPOP)-DISPM(SUBPOP):DISPM1(SUBPOP)=int(DISPM(SUBPOP)/2):DISPM2(SUBPOP)=DISPM(SUBPOP)-DISPM1(SUBPOP)
next Subpop
NSSUM=0: NMSUM=0: DispS1(NSP+1)=DispS1(1): DispM1(NSP+1)=DispM1(1): DispS2(0)=DispS2(NSP): DispM2(0)=DispM2(NSP)
for SUBPOP=1 to NSP
DRIFT=1-(1-KD/(KD+200))*rnd(1)
NS(T+1,SUBPOP)=int(NS(T+1,SUBPOP)+DISPS1(SUBPOP+1)+DISPS2(SUBPOP-1)): SNUM2=NS(T+1,SUBPOP)
NM(T+1,SUBPOP)=int(NM(T+1,SUBPOP)+DISPM1(SUBPOP+1)+DISPM2(SUBPOP-1))
if NS(T+1,SUBPOP)+NM(T+1,SUBPOP)<1 then goto *GENERATIONEND
NS(T+1,SUBPOP)=int(DRIFT*KD*NS(T+1,SUBPOP)/(NS(T+1,SUBPOP)+NM(T+1,SUBPOP)))
NM(T+1,SUBPOP)=int(DRIFT*KD*NM(T+1,SUBPOP)/(SNUM2+NM(T+1,SUBPOP)))
*GENERATIONEND
if NS(T,SUBPOP)>=1 and NS(T+1,SUBPOP)<1 then SSUBPOPEXTINCT=SSUBPOPEXTINCT+1: CUMSSUBPOPEXTINCT=CUMSSUBPOPEXTINCT+1
if NM(T,SUBPOP)>=1 and NM(T+1,SUBPOP)<1 then MSUBPOPEXTINCT=MSUBPOPEXTINCT+1: CUMMSUBPOPEXTINCT=CUMMSUBPOPEXTINCT+1
if NS(T+1,SUBPOP)+NM(T+1,SUBPOP)<1 then SubpopExtinct=SubpopExtinct+1: CUMSUBPOPEXTINCT=CUMSUBPOPEXTINCT+1
print TRIAL;" ";SUBPOP;" ";T;" ";NM(T+1,SUBPOP);" ";NS(T+1,SUBPOP);" ";MWIN;" ";SWIN;" ";: NSSUM=NSSUM+NS(T+1,SUBPOP): NMSUM=NMSUM+NM(T+1,SUBPOP)
print NMSUMatEND;" ";NSSUMatEND
next SUBPOP
NMSUMatEND=NMSUM: NSSUMatEND=NSSUM
if NSSUM=0 and NMSUM>0 then MWIN=MWIN+1: goto *TRIALEND
if NMSUM=0 and NSSUM>0 then SWIN=SWIN+1: goto *TrialEnd
if NSSUM=0 and NMSUM=0 then BOTHEXTINCT=BOTHEXTINCT+1: goto *TRIALEND
if MR=0 and SSUBPOPEXTINCT+MSUBPOPEXTINCT=NSP then goto *TRIALEND
next T
*TRIALEND
print #1,TRIAL;","; INITNS;","; INITNM;","; KD;","; MF;","; GN;","; T;","; CUMSUBPOPEXTINCT;","; CUMSSUBPOPEXTINCT;","; CUMMSUBPOPEXTINCT;","; MWIN;","; SWIN;","; BOTHEXTINCT;","; TRIAL-SWIN-MWIN-BOTHEXTINCT;","; NMSUM;",";NSSUM
SSubpopExtinct=0: MSubpopExtinct=0: TRIAL=TRIAL+1: IF TRIAL < MaxTrials+1 then GOTO *STARTTRIALS
print #1,",,,,,,,,,,poly win,mono win,population meltdown,polymorphism,,"
ENDDATE$=today$: ENDTIME=time: ENDTIME$=time$
print #1,"Time of Simulation Started= ";STARTDATE$;" ";STARTTIME$;" / Time of Simulation Ended= ";ENDDATE$;" ";ENDTIME$;" / Time for Calculation= ";ENDTIME-STARTTIME;"seconds"
close: end