Scott Reed

06/30/09

MATH 714

Traffic Circle Program

After about three hours of debugging the code for Englander’s1 traffic circle model, it no longer had compiling errors. Unfortunately, some of the subroutines it calls were not given, and so the program has linking errors that I can not resolve. The final transcript with the changes made of that program is attached.

Since that program failed to work, I began to develop my own program in FORTRAN based on the equations of Englander’straffic circle model. The transcript of the program so far is attached. So far, it completes the first two steps of the model which are determining if a car will enter the queue, if a vehicle will enter the circle, and where it will exit. The first step has been tested to show that it is working. I have begun to work on seeing if the code written does exactly what I would like it to for the second step. Based on what testing has been done, there is a deficiency in that the gap constant calculated when there are no cars in the arc preceding the entrance.

The program is titled “roundabout” and takes an input file by the same name filled with the parameters of the traffic circle, including the geometry and probabilities that are used in the program. An input file, rather than on screen inputs was chosen due to the number of these parameters, would take a considerable amount of time to enter, especially if you only want to change one or two when running a simulation. The text file is also attached, with the values given are just arbitrary at this point to develop the program. The program at this point is a bit messy, but will clean it as the program progresses.

1 Englander, Irwin. “Traffic Circle Model.” TransportationSystemsCenter. Cambridge, MA. May 1971. DOT-TSC-OST-71-6

c roundabout.f

c

c Created by: Scott Reed

c Date: June 28, 2009

c MATH 714

c

c ======

c

c Program: Roundabout

c

c Purpose: This program is used to model the flow of vehicles, on a discreet

c car level, in a traffic circle. The model is based on that

c presented by Englander*. The program reads a text file called

c "roundabout.txt" for all of the parameters of the model.

c

c*Englander, Irwin. "Traffic Circle Model." Transportation Systems

c Center. Cambridge, MA. May 1971. DOT-TSC-OST-71-6

c ======

program roundabout

c======

cVariables

c======

c NameTypeDescription

c iintegercar index

c jintegerentrance index

c lintegerexit index

c xx(j)integernumber of cars in queue at entrance j

c pxi(j)realprobability of queue increase at entrance j

c prrealrandom number between 0 and 1

c labelcharacterfile labels

c thetaj(j)realangle of intersection j in radians

c thetai(i)realangle of car i in radians

c ncariintegernumber of cars initial

c ncarfintegernumber of cars final

c kintegerprobability loop index

c nexit(i)integernumber of intersection to exit for a car i

c radiusrealradius of the traffic circle

c vel(i)realvelocity of vehicle i

c epsilonqreala constant used in the gap constant equation

c entryspdrealentry speed of a vehicle

c exitspdrealexit speed of a vehicle

c qgaprealgap constant

c thetadrealdifference in radial location of a vehicle and intersection

c thetadiffrealdifference in redial location of the closest vehicle to an intersection

c prb(k,j)realprobability of exiting at intersection k after entering at j

c sumprbrealthe sum of probabilties of exit location

implicit none

integer i,j,l,xx(4),ncari,ncarf,k,nexit(1000)

real pxi(4),pr,thetaj(5),thetai(1000),radius,vel(1000),epsilonq

real entryspd,exitspd,qgap,thetad,thetadiff,prb(4,4),sumprb

character*10 label(50)

c======

cFormat Statements

c======

1000format(a10,i3)

1100format(a10,f5.4)

1200format(a10,f9.5)

c======

cRead Circle Parameters from file

c======

open(10,file='roundabout.txt',status='old')

read(10,1200)radius

read(10,1200)epsilonq

read(10,1200)entryspd

read(10,1200)exitspd

do 100 j=1,4

read(10,1000) label(j),xx(j)

100continue

do 110 j=1,4

read(10,1100) label(j+4),pxi(j)

110continue

do 120 j=1,5

read(10,1100) label(j+8),thetaj(j)

120continue

do 130 j=1,4

do 131 k=1,4

read(10,1100) label(j+13),prb(k,j)

131 continue

130continue

print*, xx(1),xx(2),xx(3),xx(4)

print*, pxi(1),pxi(2),pxi(3),pxi(4)

print*, thetaj(1),thetaj(2),thetaj(3),thetaj(4),thetaj(5)

print*, prb(1,1),prb(2,2),prb(3,3),prb(4,4)

ncari=0

ncarf=0

c======

cstep 1: determine if a car is added to queue at intersection j

c======

do 200 j=1,4

call random_number (pr)

if (pr.gt.pxi(j)) then

xx(j)=xx(j)+1

end if

print*,pr

200continue

print*, xx(1),xx(2),xx(3),xx(4)

c======

cstep 2: determine if the first car in the queue will enter and its exit

c======

do 300 j=1,4

thetadiff=100

if (ncari.eq.0) then

go to 315

end if

do 310 i=ncari,ncarf

if (thetai(i)-thetaj(j).gt.0.0) then

if(thetai(i)-thetaj(j+1).gt.0.0) then

thetad=thetai(i)-thetaj(j)

end if

end if

if (thetad.lt.thetadiff) then

thetadiff=thetad

end if

if (thetadiff.eq.100) then

go to 315

end if

310continue

qgap= (thetadiff*radius)/(vel(i)*epsilonq)

cprint*,qgap

call random_number(pr)

if (qgap.gt.pr) then

ncari=ncari+1

vel(ncari)=entryspd

thetai(ncari)=thetaj(j)

end if

go to 305

315 ncari=ncari+1

vel(ncari)=entryspd

thetai(ncari)=thetaj(j)

go to 305

305call random_number(pr)

k=0

sumprb=0

do while (pr.lt.sumprb)

k=k+1

sumprb=sumprb+prb(k,j)

end do

nexit(ncari)=k

300continue

End

Contents of the file: roundabout.txt

radiuscir:100

epsilon_q:1.0

entrysped:5.0

exitspeed:7.0

carsinqea:2

carsinqno:2

carsinqwe:2

carsinqso:2

prbqincea:0.2

prbqincno:0.2

prbqincwe:0.2

prbqincso:0.2

thetaeast:0.0

thetanort:1.5708

thetawest:3.1416

thetasout:4.7124

thetacirc:6.2832

prbeautrn:0.25

prbearght:0.25

prbeathru:0.25

prbealeft:0.25

prbnoleft:0.25

prbnoutrn:0.25

prbnorght:0.25

prbnothru:0.25

prbwethru:0.25

prbweleft:0.25

prbweutrn:0.25

prbwerght:0.25

prbsorght:0.25

prbsothru:0.25

prbsoleft:0.25

prbsoutrn:0.25

CPROGRAM: TRAFFIC CIRCLE MODEL

CTAKEN FROM: ENGLANDER, IRWIN."TRAFFIC CIRCLE MODEL".

CTRANSPORTATIONSYSTEMSCENTER, CAMBRIDGE, MASS. 1971

CDOT-TSC-OST-71-6

CMAIN PROGRAM FOR MODEL

DIMENSION TI(51),TK(51),VI(51),VK(51),LI(51),LK(51)

DIMENSION TJ(6),XX(6),XL(6)

INTEGER QLINE,PNTR,X,Y

COMMON/A/TI,TK,VI,VK,LT,LK

COMMON/XX/X(50)

COMMON/YY/Y(50)

COMMON/PP/PNTR(2)

COMMON/NN/NCARS,NROADS

COMMON/QQ/QLINE(8)

COMMON/B/TJ,XX,XL

COMMON/C/NCAR,RC,EQ,VEN,DT

COMMON/D/ VMAX,BA,BD,ALR

COMMON/E/JX

COMMON/F/M

JX=9

PII=3.14149/180

TJ(6)=-100.

TK(51)=-100.

TI(51)=-100.

CALL SETRAN(JX)

C INPUT RC,EQ,VEN

WRITE(1,12)

11FORMAT(20H RC,VEN,EQ )

WRITE(1,11)

WRITE(1,12)

READ(1,1) RC,VEN,EQ,DT

ROTAD=RC

RC=RC/5280

13FORMAT(20H VMAX,BA,BD,ALR )

WRITE(1,13)

WRITE(1,12)

READ(1,1) VMAX,BA,BD,ALR

I=0

14FORMAT(21H THETA,VELOCITY,LEAVE)

WRITE(1,14)

2I=I+1

NCAR=I-1

WRITE(1,12)

IF(I.GT.50) GO TO 20

READ(1,9) TI(I),VI(I),LI(I)

TI(I)=TI(I)*PII

9FORMAT(2F10.5,I10)

1FORMAT(4F10.5)

12FORMAT(9X,1H.,9X,1H.,9X,1H.,9X,1H.)

IF(TI(I).GE.0.) GO TO 2

20CONTINUE

J=0

15FORMAT(20H ANGLE,QUEUE,LAMBDA )

WRITE(1,15)

4J=J+1

M=J-1

IF(J.GE.6) GO TO 1000

WRITE(1,12)

READ(1,1) TJ(J),XX(J),XL(J)

TJ(J)=TJ(J)*PII

45FORMAT(I10)

IF(TJ(J).GE.0.) GO TO 4

1000CALL ROTARY(ROTRAD,M,TJ,PNTR)

NROADS=M

10J=0

6J=J+1

IF(TJ(J).LT.0.) GO TO 7

CALL QUEUE(XX(J),XL(J))

QLINE(J)=IFIX(XX(J))

GO TO 6

7CALL ENTER

CALL VELOC

CALL POSIT

CALL ORDER

100FORMAT(6H NCARS,I5)

NCARS=NCAR

CALL CONVRT(ROTRAD,10.,TI)

CALL CHANGE

GO TO 10

END

C$1

SUBROUTINE QUEUE (X,XL)

COMMON/E/JX

CALL RAN(JX,PR)

N=X+1

XF=1.0

DO 25 I=1,N

AII=FLOAT(I)

XF=XF*AII

25CONTINUE

PX=((XL**N)*EXP(-XL))/XF

IF(PX.GE.PR) X=X+1

15FORMAT(5H QUEU,2F10.5)

RETURN

END

C$1

SUBROUTINE ENTER

DIMENSION TI(51),TK(51),VI(51),VK(51),LI(51),LK(51)

DIMENSION TJ(6),XX(6),XL(6)

COMMON/A/TI,TK,VI,VK,LI,LK

COMMON/B/TJ,XX,XL

COMMON/C/NCAR,RC,EQ,VEN,DT

COMMON/D/ VMAX,BA,BD,ALR

COMMON/E/JX

COMMON/F/M

J=1

K=1

I=1

47QJ=1

IF(TI(I).LT.0.) GO TO 17

3IF(TI(I).GT.TJ(J)) GO TO 2

15FORMAT(3I5,3F10.5)

IF(VI(I).EQ.0.) GO TO 17

QJ=((TJ(J)-TI(I))*RC)/(VI(I)*EQ)

17CALL RAN(JX,PR)

100FORMAT(3H QJ,F10.5,F10.5)

IF(QJ.LT.PR.OR.XX(J).LE.0.) GO TO 4

IF(VI(I-1).LE.0..OR.I.LE.1) GO TO 37

OJ=RC*(TI(I-1)-TJ(J))

IF(OJ.LT.0.006) GO TO 4

37TK(K)=TJ(J)

VK(K)=VEN

XX(J)=XX(J)-1.

L=0

CALL EXIT(J,L)

LK(K)=L

K=K+1

4 IF(TI(I).LT.0.) GO TO 78

CALL LEAVE(QJ,I,EQ,DT,LI(I),J)

78J=J+1

IF(J.LE.M) GO TO 47

IF(K.GT.50) RETURN

GO TO 70

2 TK(K)=TI(I)

VK(K)=VI(I)

LK(K)=LI(I)

K=K+1

I=I+1

IF(K.GE.51) RETURN

GO TO 47

70 IF(TI(I).GT.0.) GO TO 95

NCAR =K-1

TK(K)=-100.

RETURN

95TK(K)=TI(I)

LK(K)=LI(I)

VK(K)=VI(I)

K=K+1

I=I+1

GO TO 70

END

C$1

SUBROUTINE EXIT (J,L)

DIMENSION PL(5,5)

COMMON/E/JX

COMMON/F/M

AV=0.20

AB=0.20

DO 17 K=1,M

I=K+J

IF(I.GT.M)I=I-M

PL(J,K)=AV+AB

AB=AV+AB

17 CONTINUE

CALL RAN(JX,PR)

DO 1 I=1,M

L= I

IF(L.GT.M) L=L-M

IF(PR.LE.PL(J,1)) GO TO 2

1CONTINUE

2RETURN

END

C$1

SUBROUTINE LEAVE(QJ,I,EQ,DT,LII,J)

1FORMAT(4H LII,2I5)

IF(LII.NE.J) RETURN

QL=(QJ*EQ)/DT

2FORMAT(3H QL,F8.4)

IF(QL.GT.1.) RETURN

I=I+1

RETURN

END

C$1

SUBROUTINE VELOC

COMMON/C/NCAR,RC,EQ,VEN,DT

COMMON/D/ VMAX,BA,BD,ALR

COMMON/B/TJ,XX,XL

COMMON/A/TI,TK,VI,VK,LI,LK

DIMENSION TJ(6),XX(6),XL(6)

DIMENSION TI(51),TK(51),VI(51),VK(51),LI(51),LK(51)

TPI=6.28318

IF (NCAR.LE.1) GO TO 7

I=1

TP=TK(NCAR)+TPI

2 IF(TK(I).LT.0.) RETURN

LI(I)=LK(I)

VI(I)=(RC*(TP-TK(I)))/ALR

IF(VI(I).GT.VMAX) VI(I)=VMAX

DV=VI(I)-VK(I)

IF(DV.LT.BA) VI(I)=VK(I)+BA

IF(DV.GT.BD) VI(I)=VK(I)+BD

IF(VI(I).LE.0.) VI(I)=0.

TP=TK(I)

I=I+1

GO TO 2

7VI(1)=VK(1)

LI(1)=LK(1)

RETURN

END

C$1

SUBROUTINE POSIT

COMMON/C/NCAR,RC,EQ,VEN,DT

COMMON/A/TI,TK,VI,VK,LI,LK

DIMENSION TI(51),TK(51),VI(51),VK(51),LI(51),LK(51)

IF(NCAR.LE.0) RETURN

PI=6.28318

NCAR=0.

I=0

2I=I+1

IF(I.GT.50) RETURN

TI(I)=-100

IF(TK(I).LT.0.) RETURN

TI(I)=TK(I)+(VI(I)*(DT/RC))

IF(I.LE.1) GO TO 22

IF(TI(I).GT.TI(I-1)) TI(I)=TI(I-1)

22NCAR=1

1FORMAT(4H TII,F8.4)

GO TO 2

END

C$1

SUBROUTINE SETRAN (JX)

CALL RAN(JX,X)

RETURN

END

C$1

SUBROUTINE RAN(JX,YFL)

IY=JX*899

IF(IY)5,6,6

5IY=IY+32767 +1

6YFL=IY

YFL=YFL/32767.

JX=IY

RETURN

END

C$1

SUBROUTINE ORDER

COMMON/A/TI,TK,VI,VK,LI,LK

COMMON/C/NCAR,RC,EQ,VEN,DT

DIMENSION TI(51),TK(51),VI(51),VK(51),LI(51),LK(51)

IF(NCAR.LE.1) RETURN

PI=6.28318

22 IF(TI(1).GE.PI) TI(1)=TI(1)-PI

IF(TI(1).LT.0.) TI(1)=0.

IF(TI(1).GE.TI(2)) RETURN

TK(NCAR) =TI(1)

VK(NCAR)=VI(1)

LK(NCAR)=LI(1)

DO 7 J=2,NCAR

TK(J-1)=TI(J)

VK(J-1)=VT(J)

LK(J-1)=LI(J)

7CONTINUE

DO 76 J=1,NCAR

TI(J)=TK(J)

VI(J)=VK(J)

2 FORMAT(4H TIJ,F8.4)

LI(J)=LK(J)

76 CONTINUE

GO TO 22

END

C$1

SUBROUTINE ROTARY (ROTRAD,NROADS,ANGLE,PNTR)

DIMENSION ANGLE(6)

INTEGER X(6),Y(6),PNTR(2)

CDISPLAY COMMON

COMMON /DSPCOM/DSPFLG,DSPMAX/DSPUVR/DSPSZE,DSPORG,DSPINK,DSPUSR,

&DSPFLT,DSPRET,DSPTIM,DSPTOT

INTEGER DSPFLG,DSPMAX,DSPSZE,DSPORG,DSPINK,DSPUSR,DSPFLT,DSPRET,

&DSPTIM,DSPTOT

C

COMMON/DSPBFR/DSPBUF(2000)

INTEGER DSPBUF

C

CLIGHT PEN COMMON

COMMON /LPNCOM/LPNBLE,LPNBUF(50)

INTEGER LPNBLF,LPNBUF,LPNWRT

EQUIVALENCE (LPNBUF(1),LPNWRT)

C

CCHANNEL COMMON

COMMON /CHNCCM/ICHADR

C

C

CALL CLRBUF(1)

XI=ROTRAD*3.

IX=IFIX(XI)

CALL CIRCLE(512,512,IX)

DO 1 I=1,NROADS

RANG=ANGLE(1)

CRANG=COS(RANG)

SRANG=SIN(RANG)

IXDS=XI*CRANG+512.

X(I)=IXDS

IYDS=XI*SRANG+512.

Y(I)=IYDS

CALL SETPT(X(I),Y(I))

IXDS=(XI+200.)*CRANG+512.

X(I)=IXDS

IYDS=(XI+200.)*SRANG+512.

Y(I)=IYDS

CALL VECTOR(X(I),Y(I))

1CONTINUE

PNTR(1)=6+NROADS*5

DO 2 I=1,NROADS

CALL SETPT(X(I),Y(I))

CALL DSPCHR(3,2H00,2)

2 CONTINUE

PNTR(2)=PNTR(1)+NROADS*4

RETURN

END

C$1

SUBROUTINE CONVRT(ROTRAD,CARRAD,THETA)

COMMON/XX/X(50)

COMMON/YY/Y(50)

COMMON/PP/PNTR(2)

COMMON/NN/NCARS,NRAODS

DIMENSION THETA(51)

INTEGER X,Y,PNTR

2FORMAT(7h CONVRT,I5)

TRAD=(ROTRAD+CARRAD)*3.

DO 1 I=1,NCARS

RANG=THETA(I)

IDIS=TRAD*COS(RANG)+512.

X(I)=IDIS

IDIS=TRAD*SIN(RANG)+512.

1Y(I)=IDIS

RETURN

END

C$0

CRETURN

CEND

C$0

CTURN

CEND

C$0

CTURN

CEND

C$0

CRS

CRANG=THETA(I)

CIDIS=TRAD*COS(RANG)+512.

CX(I)=IDIS

CIDIS=TRAD*SIN(RANG)+512.

C 1Y(I)=IDIS

CRETURN

CEND

C&0

CRETURN

CEND

C$0

CTURN

CEND

C$0

CTURN

CEND

C$0

CRN

CEND

C$0

CTURN

CEND

C$0