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