; Vacaville Travel Demand Model
; Updated 10/15/2020 by Fehr & Peers
YEAR = 2015

;==============================================================================
;  CREATE SCENARIO NETWORK
;==============================================================================
RUN PGM=HWYNET

NETI[1]=..\Master_Network\VacaMaster.net
linki[2]=..\Input\linkcounts.csv, var=A,B,cnt_A1,cnt_P1,cnt_A2,cnt_P2,cnt_A3,cnt_P3,cnt_MD,cnt_EV,cnt_Day,cnt_DayT,Index, 
                combine=T, AVEX0=cnt_A1,cnt_P1,cnt_A2,cnt_P2,cnt_A3,cnt_P3,cnt_MD,cnt_EV,cnt_Day,cnt_DayT,Index

NETO=..\Output\?base.net, include=STNAME,DISTANCE,CAPCLASS,LANES,SPEED,cnt_A1,cnt_P1,cnt_A2,cnt_P2,cnt_A3,cnt_P3,cnt_MD,cnt_EV,cnt_Day,cnt_DayT,Index

IF (@YEAR@==2015)
  SPEED=SPEED15
  LANES=LANES15
  IF (LANES15==0)
   DELETE
  ENDIF
ENDIF

IF (@YEAR@==2035)
  SPEED=SPEED35
  LANES=LANES35
  IF (LANES35==0)
   DELETE
  ENDIF
ENDIF

IF (@YEAR@==2050)
  SPEED=SPEED50
  LANES=LANES50
  IF (LANES50==0)
   DELETE
  ENDIF
ENDIF

IF (DIST_ADJ>0) 
    DISTANCE=DIST_ADJ       ; Set gateway distance
else
    distance = sqrt((b.x-a.x)^2+(b.y-a.y)^2)/5280  ; Distance from coordinates
endif

log var=_zones

ENDRUN

;==============================================================================
;  Report out Study Intersections Correspondence
;==============================================================================
RUN PGM=HWYNET
  NETI[1]   = ..\Output\?base.net
  NODEO     = ..\Output\?_TURNS.DAT, FORMAT=TXT, INCLUDE=TURNS_OUT

  PHASE=NODEMERGE
    IF (SYNID>0|Int_ID>0)
        TURNS_OUT ='TURNS N = '+ STR(N,4,0)
    ELSE
        DELETE        
    ENDIF
  ENDPHASE
ENDRUN

;==============================================================================
;  Trip Generation
;==============================================================================
run pgm=tripgen

; Use number of zones in the highway network
zones=@hwynet._zones@

; Zonal trip generation input files

zdati[1]=..\Input\?lu.dbf,
  sum=ret,ret_hwy,service,dt,office,public,med,hosp,hosp_l,ind_h,ind_l,church,prison,elem,middle,high,college,RR,RE,AG,RetSF,RLD,RMD,RLMD,MH,RHD,RetMF

zdati[2]=..\Input\?ext.dbf, 
  sum=IX_P,XI_A,hbwp,hshp,hsrp,hscp,nhbp,
      hbwa,hsha,hsra,hsca,nhba,
      am_enter,am_exit,pm_enter,pm_exit


pao = ..\Output\?pa.DBF, DBF=YES,   ; P/A TRIP ENDS
            LIST=Z(6),P[1],P[2],P[3],P[4],P[5],P[6],P[7],P[8],P[9],P[10],P[11],P[12],P[13],P[14],P[15],P[16],P[17],P[18],P[19],P[20],
                      A[1],A[2],A[3],A[4],A[5],A[6],A[7],A[8],A[9],A[10],A[11],A[12],A[13],A[14],A[15],A[16],A[17],A[18],A[19],A[20]

LOOKUPI[1]="..\Input\triprates.CSV"
lookup lookupi=1, name=rates, 
lookup[1]=1, result=2,   ;daily
lookup[2]=1, result=3,   ;AM in
lookup[3]=1, result=4,   ;AM out
lookup[4]=1, result=5,   ;PM in
lookup[5]=1, result=6    ;PM out

; purpose split categories
array tripends=4

; combined trip categories
;     zi.1 is du, ksf, ac, students
trsf  = rates(1, 1)  *(zi.1.rr + zi.1.re + zi.1.ag + zi.1.rld + zi.1.rmd + zi.1.rlmd)
trmf  = rates(1, 2)  * zi.1.rhd   +
        rates(1, 3)  * zi.1.mh
trar  = rates(1, 4)  * zi.1.retsf +
        rates(1, 5)  * zi.1.retmf
trhwc = rates(1, 6)  * zi.1.ret_hwy
trret = rates(1, 7)  * zi.1.ret   +
        rates(1, 8)  * zi.1.service +
        rates(1, 8)  * zi.1.dt 
troff = rates(1, 9)  * zi.1.office
trpqp = rates(1,13)  * zi.1.church +
        rates(1,17)  * zi.1.public +
        rates(1,18)  * zi.1.prison +
        rates(1,22)  * zi.1.college
trmed = rates(1,14)  * zi.1.med    +
        rates(1,15)  * zi.1.hosp   +
        rates(1,16)  * zi.1.hosp_l
trind = rates(1,10)  * zi.1.ind_l  +
        rates(1,12)  * zi.1.ind_h  +
        rates(1,11)  * zi.1.hicube
trsch = rates(1,19)  * zi.1.elem   +
        rates(1,20)  * zi.1.middle +
        rates(1,21)  * zi.1.high
trhir = 0  ;not in use

loop tt=2,5
tripends[tt-1] = 
 rates(tt, 1)  *(zi.1.rr + zi.1.re + zi.1.ag + zi.1.rld + zi.1.rmd + zi.1.rlmd) +
 rates(tt, 2)  * zi.1.rhd   +
 rates(tt, 3)  * zi.1.mh    +
 rates(tt, 4)  * zi.1.retsf +
 rates(tt, 5)  * zi.1.retmf +
 rates(tt, 6)  * zi.1.ret_hwy+
 rates(tt, 7)  * zi.1.ret    +
 rates(tt, 8)  *(zi.1.dt + zi.1.service)  +
 rates(tt, 9)  * zi.1.office +
 rates(tt,10)  * zi.1.ind_l  +
 rates(tt,11)  * zi.1.hicube +
 rates(tt,12)  * zi.1.ind_h  +
 rates(tt,13)  * zi.1.church +
 rates(tt,14)  * zi.1.med    +
 rates(tt,15)  * zi.1.hosp   +
 rates(tt,16)  * zi.1.hosp_l +
 rates(tt,17)  * zi.1.public +
 rates(tt,18)  * zi.1.prison +
 rates(tt,19)  * zi.1.elem   +
 rates(tt,20)  * zi.1.middle +
 rates(tt,21)  * zi.1.high   +
 rates(tt,22)  * zi.1.college
endloop

; Lookup Tables for Production and Attraction Trip Purpose Splits
    LOOKUP, NAME=PURP,                ; Production Attraction FRACTIONS                                              
        LOOKUP[1]=1, RESULT=  2,         ; phbw
        LOOKUP[2]=1, RESULT=  3,         ; phrs
        LOOKUP[3]=1, RESULT=  4,         ; phsh
        LOOKUP[4]=1, RESULT=  5,         ; phsr
        LOOKUP[5]=1, RESULT=  6,         ; phsc
        LOOKUP[6]=1, RESULT=  7,         ; ahsrhh
        LOOKUP[7]=1, RESULT=  8,         ; pnhbhh
        LOOKUP[8]=1, RESULT=  9,         ; anhbhh
		LOOKUP[9]=1, RESULT=  10,         ; ahbw
		LOOKUP[10]=1, RESULT=  11,         ; ahrs
		LOOKUP[11]=1, RESULT=  12,         ; ahsh
		LOOKUP[12]=1, RESULT=  13,         ; ahsr
		LOOKUP[13]=1, RESULT=  14,         ; ahsc
		LOOKUP[14]=1, RESULT=  15,         ; pnhb
		LOOKUP[15]=1, RESULT=  16,         ; anhb
        INTERPOLATE=FALSE,                     
        FILE="..\Input\PA_splits.CSV"

; purpose splits (daily)
phbw	= PURP(1, 1) * trsf  +  PURP(1, 2) * trmf  +  PURP(1, 3) * trar
phrs	= PURP(2, 1) * trsf  +  PURP(2, 2) * trmf  +  PURP(2, 3) * trar  ;regional shopping (hwc, outlets)
phsh	= PURP(3, 1) * trsf  +  PURP(3, 2) * trmf  +  PURP(3, 3) * trar  ;local shopping
phsr	= PURP(4, 1) * trsf  +  PURP(4, 2) * trmf  +  PURP(4, 3) * trar
phsc	= PURP(5, 1) * trsf  +  PURP(5, 2) * trmf  +  PURP(5, 3) * trar
ahsrhh	= PURP(6, 1) * trsf  +  PURP(6, 2) * trmf  +  PURP(6, 3) * trar
pnhbhh	= PURP(7, 1) * trsf  +  PURP(7, 2) * trmf  +  PURP(7, 3) * trar
anhbhh	= PURP(8, 1) * trsf  +  PURP(8, 2) * trmf  +  PURP(8, 3) * trar

ahbw	= PURP(9, 4) * trhwc  +  PURP(9, 5) * trret  +  PURP(9, 6) * troff  +  PURP(9, 7) * trpqp  +  PURP(9, 8) * trmed  +  PURP(9, 9) * trind  +  PURP(9, 10) * trsch  +  PURP(9, 11) * trhir
ahrs	= PURP(10, 4) * trhwc  +  PURP(10, 5) * trret  +  PURP(10, 6) * troff  +  PURP(10, 7) * trpqp  +  PURP(10, 8) * trmed  +  PURP(10, 9) * trind  +  PURP(10, 10) * trsch  +  PURP(10, 11) * trhir
ahsh	= PURP(11, 4) * trhwc  +  PURP(11, 5) * trret  +  PURP(11, 6) * troff  +  PURP(11, 7) * trpqp  +  PURP(11, 8) * trmed  +  PURP(11, 9) * trind  +  PURP(11, 10) * trsch  +  PURP(11, 11) * trhir
ahsr	= PURP(12, 4) * trhwc  +  PURP(12, 5) * trret  +  PURP(12, 6) * troff  +  PURP(12, 7) * trpqp  +  PURP(12, 8) * trmed  +  PURP(12, 9) * trind  +  PURP(12, 10) * trsch  +  PURP(12, 11) * trhir
ahsc	= PURP(13, 4) * trhwc  +  PURP(13, 5) * trret  +  PURP(13, 6) * troff  +  PURP(13, 7) * trpqp  +  PURP(13, 8) * trmed  +  PURP(13, 9) * trind  +  PURP(13, 10) * trsch  +  PURP(13, 11) * trhir
pnhb	= PURP(14, 4) * trhwc  +  PURP(14, 5) * trret  +  PURP(14, 6) * troff  +  PURP(14, 7) * trpqp  +  PURP(14, 8) * trmed  +  PURP(14, 9) * trind  +  PURP(14, 10) * trsch  +  PURP(14, 11) * trhir
anhb	= PURP(15, 4) * trhwc  +  PURP(15, 5) * trret  +  PURP(15, 6) * troff  +  PURP(15, 7) * trpqp  +  PURP(15, 8) * trmed  +  PURP(15, 9) * trind  +  PURP(15, 10) * trsch  +  PURP(15, 11) * trhir


; Set Up Lookup Tables for IX and XI fractions
    LOOKUP, NAME=IX,                     ; IX FRACTIONS                                                   
        LOOKUP[1]=1, RESULT=  2,         ; HBW               
        LOOKUP[2]=1, RESULT=  3,         ; HRS               
        LOOKUP[3]=1, RESULT=  4,         ; HSH               
        LOOKUP[4]=1, RESULT=  5,         ; HSR            
        LOOKUP[5]=1, RESULT=  6,         ; HSC              
        LOOKUP[6]=1, RESULT=  7,         ; NBH               
        LOOKUP[7]=1, RESULT=  14,         ; AM              
        LOOKUP[8]=1, RESULT=  15,         ; PM  
        INTERPOLATE=FALSE,                     
        FILE="..\Input\IXXI_Percentages.CSV"
	
	LOOKUP, NAME=XI,                     ; IX FRACTIONS                                                   
        LOOKUP[1]=1, RESULT=  8,         ; HBW               
        LOOKUP[2]=1, RESULT=  9,         ; HRS               
        LOOKUP[3]=1, RESULT=  10,         ; HSH               
        LOOKUP[4]=1, RESULT=  11,         ; HSR            
        LOOKUP[5]=1, RESULT=  12,         ; HSC              
        LOOKUP[6]=1, RESULT=  13,         ; NBH            
        LOOKUP[7]=1, RESULT=  16,         ; AM              
        LOOKUP[8]=1, RESULT=  17,         ; PM  		
        INTERPOLATE=FALSE,                     
        FILE="..\Input\IXXI_Percentages.CSV"
                                                          

; Raw vehicle trip generation (daily, P-A by purpose)
; II PAs
p[1] = (phbw*(1-IX(1, I)))
p[2] = (phrs*(1-IX(2, I)))
p[3] = (phsh*(1-IX(3, I)))
p[4] = (phsr*(1-IX(4, I)))
p[5] = (phsc*(1-IX(5, I)))
p[6] = ((pnhb + pnhbhh)*(1-IX(6, I)))

a[1] = (ahbw*(1-XI(1, I)))
a[2] = (ahrs*(1-XI(2, I)))
a[3] = (ahsh*(1-XI(3, I)))
a[4] = ((ahsr + ahsrhh)*(1-XI(4, I)))
a[5] = (ahsc*(1-XI(5, I)))
a[6] = ((anhb + anhbhh)*(1-XI(6, I)))


; IX PAs
; Raw vehicle trip generation (daily, P-A by purpose)
p[9]  = (phbw*IX(1, I))   
p[11] = (phrs*IX(2, I))      
p[13] = (phsh*IX(3, I))         
p[15] = (phsr*IX(4, I))          
p[17] = (phsc*IX(5, I))         
p[19] = ((pnhb + pnhbhh)*IX(6, I)) 

a[9]  = zi.2.XI_A
a[11] = zi.2.XI_A
a[13] = zi.2.XI_A
a[15] = zi.2.XI_A
a[17] = zi.2.XI_A
a[19] = zi.2.XI_A

; XI PAs
a[10] = (ahbw*XI(1, I))   
a[12] = (ahrs*XI(2, I))     
a[14] = (ahsh*XI(3, I))       
a[16] = ((ahsr + ahsrhh)*XI(4, I))
a[18] = (phsc*XI(5, I))	 
a[20] = ((anhb + anhbhh)*XI(6, I)) 

P[10] = zi.2.IX_P 
P[12] = zi.2.IX_P
P[14] = zi.2.IX_P
P[16] = zi.2.IX_P 
P[18] = zi.2.IX_P 
P[20] = zi.2.IX_P 


; Peak OD vehicle trip generation
p[7]  = tripends[2] + tripends[2]*(IX(7, I))
p[8]  = tripends[4] + tripends[4]*(IX(8, I))
a[7]  = tripends[1] + tripends[1]*(XI(7, I))
a[8]  = tripends[3] + tripends[3]*(XI(8, I))


phase=adjust

; Print PAs before Balancing
PTOT1_U = PTOT(1)
PTOT2_U = PTOT(2)
PTOT3_U = PTOT(3)
PTOT4_U = PTOT(4)
PTOT5_U = PTOT(5)
PTOT6_U = PTOT(6)
PTOT7_U = PTOT(7)
PTOT8_U = PTOT(8)
PTOT9_U = PTOT(9)
PTOT10_U = PTOT(10)
PTOT11_U = PTOT(11)
PTOT12_U = PTOT(12)
PTOT13_U = PTOT(13)
PTOT14_U = PTOT(14)
PTOT15_U = PTOT(15)
PTOT16_U = PTOT(16)
PTOT17_U = PTOT(17)
PTOT18_U = PTOT(18)
PTOT19_U = PTOT(19)
PTOT20_U = PTOT(20)
;PTOT21_U = PTOT(21)
;PTOT22_U = PTOT(22)

ATOT1_U = ATOT(1)
ATOT2_U = ATOT(2)
ATOT3_U = ATOT(3)
ATOT4_U = ATOT(4)
ATOT5_U = ATOT(5)
ATOT6_U = ATOT(6)
ATOT7_U = ATOT(7)
ATOT8_U = ATOT(8)
ATOT9_U = ATOT(9)
ATOT10_U = ATOT(10)
ATOT11_U = ATOT(11)
ATOT12_U = ATOT(12)
ATOT13_U = ATOT(13)
ATOT14_U = ATOT(14)
ATOT15_U = ATOT(15)
ATOT16_U = ATOT(16)
ATOT17_U = ATOT(17)
ATOT18_U = ATOT(18)
ATOT19_U = ATOT(19)
ATOT20_U = ATOT(20)
;ATOT21_U = ATOT(21)
;ATOT22_U = ATOT(22)

; ***P/A Balancing***

; II PAs
a[1] = a[1] / a[1][0] * p[1][0]
a[2] = a[2] / a[2][0] * p[2][0]
a[3] = a[3] / a[3][0] * p[3][0]
a[4] = a[4] / a[4][0] * p[4][0]
a[5] = a[5] / a[5][0] * p[5][0]
p[6] = (p[6] + a[6])/2
a[6] = (p[6] + a[6])/2
p[7] = p[7] / p[7][0] * (p[7][0] + a[7][0]) / 2
a[7] = a[7] / a[7][0] * (p[7][0] + a[7][0]) / 2
p[8] = p[8] / p[8][0] * (p[8][0] + a[8][0]) / 2
a[8] = a[8] / a[8][0] * (p[8][0] + a[8][0]) / 2

;IX/XI PAs
a[9]  = a[9]  / a[9][0] * p[9][0]
a[11] = a[11] / a[11][0] * p[11][0]
a[13] = a[13] / a[13][0] * p[13][0]
a[15] = a[15] / a[15][0] * p[15][0]
a[17] = a[17] / a[17][0] * p[17][0]
A[19] = (A[19] + P[19])/2

P[10] = P[10] / P[10][0] * A[10][0]
P[12] = P[12] / P[12][0] * A[12][0]
P[14] = P[14] / P[14][0] * A[14][0]
P[16] = P[16] / P[16][0] * A[16][0]
P[18] = P[18] / P[18][0] * A[18][0]
P[20] = (A[20] + P[20])/2

;p[21] = p[21] / p[21][0] * (p[21][0] + a[21][0]) / 2
;a[21] = a[21] / a[21][0] * (p[21][0] + a[21][0]) / 2
;p[22] = p[22] / p[22][0] * (p[22][0] + a[22][0]) / 2
;a[22] = a[22] / a[22][0] * (p[22][0] + a[22][0]) / 2

; Print PAs After Balancing
PTOT1_B = PTOT(1)
PTOT2_B = PTOT(2)
PTOT3_B = PTOT(3)
PTOT4_B = PTOT(4)
PTOT5_B = PTOT(5)
PTOT6_B = PTOT(6)
PTOT7_B = PTOT(7)
PTOT8_B = PTOT(8)
PTOT9_B = PTOT(9)
PTOT10_B = PTOT(10)
PTOT11_B = PTOT(11)
PTOT12_B = PTOT(12)
PTOT13_B = PTOT(13)
PTOT14_B = PTOT(14)
PTOT15_B = PTOT(15)
PTOT16_B = PTOT(16)
PTOT17_B = PTOT(17)
PTOT18_B = PTOT(18)
PTOT19_B = PTOT(19)
PTOT20_B = PTOT(20)
;PTOT21_B = PTOT(21)
;PTOT22_B = PTOT(22)

ATOT1_B = ATOT(1)
ATOT2_B = ATOT(2)
ATOT3_B = ATOT(3)
ATOT4_B = ATOT(4)
ATOT5_B = ATOT(5)
ATOT6_B = ATOT(6)
ATOT7_B = ATOT(7)
ATOT8_B = ATOT(8)
ATOT9_B = ATOT(9)
ATOT10_B = ATOT(10)
ATOT11_B = ATOT(11)
ATOT12_B = ATOT(12)
ATOT13_B = ATOT(13)
ATOT14_B = ATOT(14)
ATOT15_B = ATOT(15)
ATOT16_B = ATOT(16)
ATOT17_B = ATOT(17)
ATOT18_B = ATOT(18)
ATOT19_B = ATOT(19)
ATOT20_B = ATOT(20)
;ATOT21_B = ATOT(21)
;ATOT22_B = ATOT(22)

PRINT FORM=13,      
    LIST="BEFORE BALANCING","\n",                      
         ",,HBW,HRS,HSH,HSR,HSC,NHB,AM,PM,HBW_IX, HBW_XI, HRS_IX, HRS_XI, HSH_IX, HSH_XI, HSR_I, HSR_I, HSC_IX, HSC_XI, NHB_IX, NHB_XI","\n",                      
         ",PTOT,",PTOT1_U,",",PTOT2_U,",",PTOT3_U,",",PTOT4_U,",",PTOT5_U,",",PTOT6_U,",",PTOT7_U,",",PTOT8_U,",",PTOT9_U,",",PTOT10_U,",",PTOT11_U,",",PTOT12_U,",",PTOT13_U,",",PTOT14_U,",",PTOT15_U,",",PTOT16_U,",",PTOT17_U,",",PTOT18_U,",",PTOT19_U,",",PTOT20_U,"\n",                      
         ",ATOT,",ATOT1_U,",",ATOT2_U,",",ATOT3_U,",",ATOT4_U,",",ATOT5_U,",",ATOT6_U,",",ATOT7_U,",",ATOT8_U,",",ATOT9_U,",",ATOT10_U,",",ATOT11_U,",",ATOT12_U,",",ATOT13_U,",",ATOT14_U,",",ATOT15_U,",",ATOT16_U,",",ATOT17_U,",",ATOT18_U,",",ATOT19_U,",",ATOT20_U,"\n",                      
		"AFTER BALANCING","\n",                      
         ",,HBW,HRS,HSH,HSR,HSC,NHB,AM,PM,HBW_IX, HBW_XI, HRS_IX, HRS_XI, HSH_IX, HSH_XI, HSR_I, HSR_I, HSC_IX, HSC_XI, NHB_IX, NHB_XI","\n",                      
         ",PTOT,",PTOT1_B,",",PTOT2_B,",",PTOT3_B,",",PTOT4_B,",",PTOT5_B,",",PTOT6_B,",",PTOT7_B,",",PTOT8_B,",",PTOT9_B,",",PTOT10_B,",",PTOT11_B,",",PTOT12_B,",",PTOT13_B,",",PTOT14_B,",",PTOT15_B,",",PTOT16_B,",",PTOT17_B,",",PTOT18_B,",",PTOT19_B,",",PTOT20_B,"\n",                      
         ",ATOT,",ATOT1_B,",",ATOT2_B,",",ATOT3_B,",",ATOT4_B,",",ATOT5_B,",",ATOT6_B,",",ATOT7_B,",",ATOT8_B,",",ATOT9_B,",",ATOT10_B,",",ATOT11_B,",",ATOT12_B,",",ATOT13_B,",",ATOT14_B,",",ATOT15_B,",",ATOT16_B,",",ATOT17_B,",",ATOT18_B,",",ATOT19_B,",",ATOT20_B,"\n",    
    FILE="..\Output\?_PA_Summary.CSV", REWIND=T  

ENDRUN

;==============================================================================
;  CREATE SKIM MATRIX
;==============================================================================

RUN PGM=HWYLOAD

NETI=..\Output\?base.net
FILEI TURNPENI =..\Input\TURNPEN.DAT, MISSINGLINK=1  ; INCLUDE TURN PENALTIES
MATO=..\Output\Skim.mat, MO=1, Name=IMP

PHASE=ILOOP
  Path=TIME, PENI=1,
  MW[1]=PATHTRACE(TIME,1), noaccess=0
  MW[1][I]=LOWEST(1,1,0.001,999,I)/2
; No access for XX movements
  if (i=1-11)
      jloop j=1,11
          mw[1]=0
      endjloop
  endif

  ;Terminal Times
  jloop
      if (mw[1] > 0) mw[1] = mw[1] + 2  ;constant 1 min each trip end
  endjloop
ENDPHASE

ENDRUN

;==============================================================================
;  TRIP DISTRIBUTION
;==============================================================================

RUN PGM=TRIPDIST

FILEI MATI[1] = ..\Output\Skim.MAT
;  FILEI MATI[2] = ..\Input\KFAC.MAT
FILEI LOOKUPI[1] = ..\Input\Vacaville_FF.csv
FILEI ZDATI[1] = ..\Output\?pa.DBF, Z=Z 

;'Z=#1, Zone=#1, P1=#2, P2=#3, P3=#4, P4=#5, P5=#6, P6=#7, P7=#8, P8=#9, P9=#10, P10=#11, P11=#12, P12=#13, P13=#14, P14=#15, P15=#16, P16=#17, A1=#18, A2=#19, A3=#20, A4=#21, A5=#22, A6=#23, A7=#24, A8=#25, A9=#24, A10=#25, A11=#26, A12=#27, A13=#28, A14=#29, A15=#30, A16=#31'
                                        
                                         
FILEO MATO[1] = ..\Output\TripDist.MAT, MO=1-19, name=HBW,HRS,HSH,HSR,HSC,NHB,HBW_IX, HRS_IX, HSH_IX, HSR_I, HSC_IX, NHB_IX, HBW_XI, HRS_XI, HSH_XI, HSR_I, HSC_XI, NHB_XI,TOTAL                ;Daily P-A trip matrix

  MAXITERS=99
  MAXRMSE = 0.01
  ZONEMSG= 100


;FRICTION FACTOR LOOKUP
  LOOKUP LOOKUPI=1, NAME=FF,
    LOOKUP[1]=1,RESULT=2,                 ;HBW
    LOOKUP[2]=1,RESULT=4,                 ;HRS
    LOOKUP[3]=1,RESULT=3,                 ;HSH
    LOOKUP[4]=1,RESULT=4,                 ;HSR
    LOOKUP[5]=1,RESULT=5,                 ;HSC
    LOOKUP[6]=1,RESULT=6,                 ;NHB
	LOOKUP[7]=1,RESULT=2,                 ;HBW
    LOOKUP[8]=1,RESULT=4,                 ;HRS
    LOOKUP[9]=1,RESULT=3,                 ;HSH
    LOOKUP[10]=1,RESULT=4,                 ;HSR
    LOOKUP[11]=1,RESULT=5,                 ;HSC
    LOOKUP[12]=1,RESULT=6,                 ;NHB
	LOOKUP[13]=1,RESULT=2,                 ;HBW
    LOOKUP[14]=1,RESULT=4,                 ;HRS
    LOOKUP[15]=1,RESULT=3,                 ;HSH
    LOOKUP[16]=1,RESULT=4,                 ;HSR
    LOOKUP[17]=1,RESULT=5,                 ;HSC
    LOOKUP[18]=1,RESULT=6,                 ;NHB
	INTERPOLATE=Y, SETUPPER=N

;SETUP THE WORKING P AND A
  SETPA P[1]=P1  	A[1]=A1
  SETPA P[2]=P2  	A[2]=A2
  SETPA P[3]=P3  	A[3]=A3
  SETPA P[4]=P4  	A[4]=A4
  SETPA P[5]=P5  	A[5]=A5
  SETPA P[6]=P6  	A[6]=A6
  SETPA P[7]=P9  	A[7]=A9
  SETPA P[8]=P10  	A[8]=A10
  SETPA P[9]=P11	A[9]=A11
  SETPA P[10]=P12  	A[10]=A12
  SETPA P[11]=P13  	A[11]=A13
  SETPA P[12]=P14  	A[12]=A14
  SETPA P[13]=P15  	A[13]=A15
  SETPA P[14]=P16  	A[14]=A16
  SETPA P[15]=P17  	A[15]=A17
  SETPA P[16]=P18  	A[16]=A18
  SETPA P[17]=P19  	A[17]=A19
  SETPA P[18]=P20  	A[18]=A20
  
  mw[21] = mi.1.1
  if (i=1-11)
      jloop j=1,11
         mw[21] = -100 							; WHAT DOES THE -100 DO?
      endjloop
  endif

  GRAVITY   PURPOSE=1, LOS=mw[21], FFACTORS=FF  ; HBW
  GRAVITY   PURPOSE=2, LOS=mw[21], FFACTORS=FF  ; HRS
  GRAVITY   PURPOSE=3, LOS=mw[21], FFACTORS=FF  ; HSH
  GRAVITY   PURPOSE=4, LOS=mw[21], FFACTORS=FF  ; HSR
  GRAVITY   PURPOSE=5, LOS=mw[21], FFACTORS=FF  ; HSC
  GRAVITY   PURPOSE=6, LOS=mw[21], FFACTORS=FF  ; NHB
  GRAVITY   PURPOSE=7, LOS=mw[21], FFACTORS=FF  ; HBW
  GRAVITY   PURPOSE=8, LOS=mw[21], FFACTORS=FF  ; HRS
  GRAVITY   PURPOSE=9, LOS=mw[21], FFACTORS=FF  ; HSH
  GRAVITY   PURPOSE=10, LOS=mw[21], FFACTORS=FF  ; HSR
  GRAVITY   PURPOSE=11, LOS=mw[21], FFACTORS=FF  ; HSC
  GRAVITY   PURPOSE=12, LOS=mw[21], FFACTORS=FF  ; NHB
  GRAVITY   PURPOSE=13, LOS=mw[21], FFACTORS=FF  ; HBW
  GRAVITY   PURPOSE=14, LOS=mw[21], FFACTORS=FF  ; HRS
  GRAVITY   PURPOSE=15, LOS=mw[21], FFACTORS=FF  ; HSH
  GRAVITY   PURPOSE=16, LOS=mw[21], FFACTORS=FF  ; HSR
  GRAVITY   PURPOSE=17, LOS=mw[21], FFACTORS=FF  ; HSC
  GRAVITY   PURPOSE=18, LOS=mw[21], FFACTORS=FF  ; NHB
  
  MW[19]=MW[1]+MW[2]+MW[3]+MW[4]+MW[5]+MW[6]+MW[7]+MW[8]+MW[9]+MW[10]+MW[11]+MW[12]+MW[13]+MW[14]+MW[15]+MW[16]+MW[17]+MW[18]  ;Total

FREQUENCY BASEMW=21,VALUEMW=1,RANGE=1-50, TITLE='Work Trips'
FREQUENCY BASEMW=21,VALUEMW=2,RANGE=1-50, TITLE='Regional & Hwy Shop'
FREQUENCY BASEMW=21,VALUEMW=3,RANGE=1-50, TITLE='Local Shop'
FREQUENCY BASEMW=21,VALUEMW=4,RANGE=1-50, TITLE='Social-Recreational'
FREQUENCY BASEMW=21,VALUEMW=5,RANGE=1-50, TITLE='School'
FREQUENCY BASEMW=21,VALUEMW=6,RANGE=1-50, TITLE='Non Home-Based'

FREQUENCY BASEMW=21,VALUEMW=7,RANGE=1-50, TITLE='Work Trips IX'
FREQUENCY BASEMW=21,VALUEMW=8,RANGE=1-50, TITLE='Regional & Hwy Shop IX'
FREQUENCY BASEMW=21,VALUEMW=9,RANGE=1-50, TITLE='Local Shop IX'
FREQUENCY BASEMW=21,VALUEMW=10,RANGE=1-50, TITLE='Social-Recreational IX'
FREQUENCY BASEMW=21,VALUEMW=11,RANGE=1-50, TITLE='School IX'
FREQUENCY BASEMW=21,VALUEMW=12,RANGE=1-50, TITLE='Non Home-Based IX'

FREQUENCY BASEMW=21,VALUEMW=13,RANGE=1-50, TITLE='Work Trips XI'
FREQUENCY BASEMW=21,VALUEMW=14,RANGE=1-50, TITLE='Regional & Hwy Shop XI'
FREQUENCY BASEMW=21,VALUEMW=15,RANGE=1-50, TITLE='Local Shop XI'
FREQUENCY BASEMW=21,VALUEMW=16,RANGE=1-50, TITLE='Social-Recreational XI'
FREQUENCY BASEMW=21,VALUEMW=17,RANGE=1-50, TITLE='School XI'
FREQUENCY BASEMW=21,VALUEMW=18,RANGE=1-50, TITLE='Non Home-Based XI'

ENDRUN


;==============================================================================
;  XX Trips
;==============================================================================

RUN PGM=MATRIX

MATI=..\Input\?XXda.csv, pattern=IJ:V FIELDS=#1,0,2-12
ZONES=11
MATO=..\Output\?XXda.mat, MO=1

MW[1]=MI.1.1

ENDRUN

;==============================================================================
;  Basic P-A to peak period O-D conversion
;==============================================================================

RUN PGM=MATRIX 

    MATI[1] = ..\Output\TripDist.mat     ;  PA TRIP MATRICES
    mati[2] = ..\Output\?xxda.mat         ; thru trips

    MATO[1] = ..\Output\A3_OD_ByPurp.MAT, mo=21-28,dec=8*4, name=HBW,HRS,HSH,HSR,HSC,NHB,XX,Total   ; A3 peak period trip table
    MATO[2] = ..\Output\MD_OD_ByPurp.MAT, mo=31-38,dec=8*4, name=HBW,HRS,HSH,HSR,HSC,NHB,XX,Total   ; MD peak period trip table
    MATO[3] = ..\Output\P3_OD_ByPurp.MAT, mo=41-48,dec=8*4, name=HBW,HRS,HSH,HSR,HSC,NHB,XX,Total   ; P3 peak period trip table
    MATO[4] = ..\Output\EV_OD_ByPurp.MAT, mo=51-58,dec=8*4, name=HBW,HRS,HSH,HSR,HSC,NHB,XX,Total   ; EV peak period trip table
    MATO[5] = ..\Output\A1_OD_Temp.MAT,   mo=61-68,dec=8*4, name=HBW,HRS,HSH,HSR,HSC,NHB,XX,Total   ; A1 peak period trip table
    MATO[6] = ..\Output\P1_OD_Temp.MAT,   mo=71-78,dec=8*4, name=HBW,HRS,HSH,HSR,HSC,NHB,XX,Total   ; P1 peak period trip table
;

    LOOKUP,                             ; PEAK FACTOR LOOKUP
        FILE="..\Input\todfactors.CSV",               
        NAME=PHFAC,                      
        LOOKUP[1] =1, RESULT= 2,        ; DEP_HBW
        LOOKUP[2] =1, RESULT= 3,        ; RET_HBW
        LOOKUP[3] =1, RESULT= 4,        ; DEP_HRS
        LOOKUP[4] =1, RESULT= 5,        ; RET_HRS
        LOOKUP[5] =1, RESULT= 4,        ; DEP_HSH
        LOOKUP[6] =1, RESULT= 5,        ; RET_HSH
        LOOKUP[7] =1, RESULT= 6,        ; DEP_HSR
        LOOKUP[8] =1, RESULT= 7,        ; RET_HSR
        LOOKUP[9] =1, RESULT= 8,        ; DEP_HSC
        LOOKUP[10]=1, RESULT= 9,        ; RET_HSC
        LOOKUP[11]=1, RESULT=10,        ; DEP_NHB
        LOOKUP[12]=1, RESULT=11,        ; RET_NHB
        LOOKUP[13]=1, RESULT=12,        ; DEP_XX
        LOOKUP[14]=1, RESULT=13         ; RET_XX

;-------------------------------------------------------------------------------
    ;
    ;  PERIOD TRIPS
    
;  p to a direction   (number internal matrices same as factor-lookup index)
    mw[01] = mi.1.1 + mi.1.7  + mi.1.13
	mw[03] = mi.1.2 + mi.1.8  + mi.1.14
    mw[05] = mi.1.3 + mi.1.9  + mi.1.15
    mw[07] = mi.1.4 + mi.1.10 + mi.1.16
    mw[09] = mi.1.5 + mi.1.11 + mi.1.17
    mw[11] = mi.1.6 + mi.1.12 + mi.1.18
    mw[13] = mi.2.1

;  a to p direction
    mw[02] = mi.1.1.t + mi.1.7.t  + mi.1.13.t
    mw[04] = mi.1.2.t + mi.1.8.t  + mi.1.14.t
    mw[06] = mi.1.3.t + mi.1.9.t  + mi.1.15.t
    mw[08] = mi.1.4.t + mi.1.10.t + mi.1.16.t
    mw[10] = mi.1.5.t + mi.1.11.t + mi.1.17.t
    mw[12] = mi.1.6.t + mi.1.12.t + mi.1.18.t
    mw[14] = mi.2.1.t
   
;   convert lookups to array for efficiency
    array a3facs=16, mdfacs=16, p3facs=16, evfacs=16, a1facs=16, p1facs=16
    if (i=1)
        loop jj=1,14
            a3facs[jj] = PHFAC(jj,1)
            mdfacs[jj] = PHFAC(jj,2)
            p3facs[jj] = PHFAC(jj,3)
            evfacs[jj] = PHFAC(jj,4)
            a1facs[jj] = PHFAC(jj,5)
            p1facs[jj] = PHFAC(jj,6)
        endloop
    endif
;-------------------------------------------------------------------------------

    loop jj=1,zones   ;(for loop inside jloop)
        a3trips = 0
        mdtrips = 0
        p3trips = 0
        evtrips = 0
        a1trips = 0
        p1trips = 0
        loop purp=1,7
            luap = purp + purp   ; 2,4,6...14
            lupa = luap - 1      ; 1,3,5...13
            a3mw = purp + 20
            mdmw = purp + 30
            p3mw = purp + 40
            evmw = purp + 50
            a1mw = purp + 60
            p1mw = purp + 70
            a3purptrips = mw[lupa][jj] * a3facs[lupa] + mw[luap][jj] * a3facs[luap]
            mdpurptrips = mw[lupa][jj] * mdfacs[lupa] + mw[luap][jj] * mdfacs[luap]
            p3purptrips = mw[lupa][jj] * p3facs[lupa] + mw[luap][jj] * p3facs[luap]
            evpurptrips = mw[lupa][jj] * evfacs[lupa] + mw[luap][jj] * evfacs[luap]
            a1purptrips = mw[lupa][jj] * a1facs[lupa] + mw[luap][jj] * a1facs[luap]
            p1purptrips = mw[lupa][jj] * p1facs[lupa] + mw[luap][jj] * p1facs[luap]
            mw[a3mw][jj] = a3purptrips
            mw[mdmw][jj] = mdpurptrips
            mw[p3mw][jj] = p3purptrips
            mw[evmw][jj] = evpurptrips
            mw[a1mw][jj] = a1purptrips
            mw[p1mw][jj] = p1purptrips
            a3trips = a3trips + a3purptrips
            mdtrips = mdtrips + mdpurptrips
            p3trips = p3trips + p3purptrips
            evtrips = evtrips + evpurptrips
            a1trips = a1trips + a1purptrips
            p1trips = p1trips + p1purptrips
        endloop
        mw[28][jj] = a3trips
        mw[38][jj] = mdtrips
        mw[48][jj] = p3trips
        mw[58][jj] = evtrips
        mw[68][jj] = a1trips
        mw[78][jj] = p1trips
    endloop
 
ENDRUN

;==========================================================================
; Match peak period trip ends
;==========================================================================

RUN PGM=FRATAR ;MSG='Factor vehicle trips to match original trips in and out of TAZs'

FILEI MATI[1] = ..\Output\A1_OD_Temp.MAT
FILEI MATI[2] = ..\Output\P1_OD_Temp.MAT
FILEI zdati[1] = ..\Output\?PA.dbf, z=z
;'z=#1, AMO=#8, PMO=#9, AMOXI=#17, PMOIX=#18, AMD=#24, PMD=#25, AMDXI=#30, PMDIX=#31'

MATO[1]= ..\Output\A1_OD_Tot.MAT, MO=1, dec=4, NAME=AM_Total
MATO[2]= ..\Output\P1_OD_Tot.MAT, MO=2, dec=4, NAME=PM_Total

; Factor trip tables to match original ins and outs from unadjusted trip tables
MAXRMSE=0.01 MAXITERS=10

SETPA P[1]=P7, A[1]=A7, MW[1]=MI.1.Total
SETPA P[2]=P8, A[2]=A8, MW[2]=MI.2.Total 

ENDRUN

;==========================================================================
; Split peak period trips by purpose
;==========================================================================
run pgm=matrix
mati[1] = ..\Output\A1_OD_Temp.mat
mati[2] = ..\Output\P1_OD_Temp.mat

mati[3] = ..\Output\A1_OD_Tot.mat
mati[4] = ..\Output\P1_OD_Tot.mat

mato[1] = ..\Output\A1_OD_ByPurp.mat, mo=1-8,   dec=8*3, name=HBW,HRS,HSH,HSR,HSC,NHB,XX,Total
mato[2] = ..\Output\P1_OD_ByPurp.mat, mo=11-18, dec=8*3, name=HBW,HRS,HSH,HSR,HSC,NHB,XX,Total

jloop
   if (mi.1.Total > 0)
      ratio = mi.3.AM_Total / mi.1.Total
      mw[01] = mi.1.1 * ratio
      mw[02] = mi.1.2 * ratio
      mw[03] = mi.1.3 * ratio
      mw[04] = mi.1.4 * ratio
      mw[05] = mi.1.5 * ratio
      mw[06] = mi.1.6 * ratio
      mw[07] = mi.1.7 * ratio
      mw[08] = mi.3.AM_Total
   endif
   if (mi.2.Total > 0)
      ratio = mi.4.PM_Total / mi.2.Total
      mw[11] = mi.2.1 * ratio
      mw[12] = mi.2.2 * ratio
      mw[13] = mi.2.3 * ratio
      mw[14] = mi.2.4 * ratio
      mw[15] = mi.2.5 * ratio
      mw[16] = mi.2.6 * ratio
      mw[17] = mi.2.7 * ratio
      mw[18] = mi.4.PM_Total
   endif
endjloop

endrun


;==============================================================================
;  AM PEAK PERIOD ASSIGNMENT
;==============================================================================
;==============================================================================
RUN PGM=HWYLOAD

NETI=..\Output\?base.net
FILEI TURNPENI = ..\Input\TURNPEN.DAT, MISSINGLINK=1
FILEI MATI[1] = ..\Output\A3_OD_ByPurp.MAT
FILEO NETO = ..\Output\?_A3.net

TURNVOLO[1]=..\Output\?_A3_turn.dbf,FORMAT=DBF, DEC=0

;-----------------------------------------------------------------------------;
  ; CAPACITY IN VPH/LANE:
  SPDCAP CAPACITY[1]=2000       ; FREEWAY
  SPDCAP CAPACITY[2]=1000       ; RURAL ROAD
  SPDCAP CAPACITY[3]=800        ; MAJOR ARTERIAL
  SPDCAP CAPACITY[4]=700        ; MINOR ARTERIAL
  SPDCAP CAPACITY[5]=600        ; COLLECTOR
  SPDCAP CAPACITY[6]=1500       ; LINEAR RAMP
  SPDCAP CAPACITY[7]=1200       ; LOOP RAMP
  SPDCAP CAPACITY[10]=9999      ; CENTROID CONNECTOR
;  SPDCAP CAPACITY[10]=3000      ; CENTROID CONNECTOR
;-----------------------------------------------------------------------------;
  PARAMETERS CAPFAC=2.85  ;3-hour intensity factor
  PARAMETERS MAXITERS=70
  parameters relativegap=0.0002
  RMSE=0
  PARAMETERS GAP = 0
  PARAMETERS RAAD = 0
  PARAMETERS AAD = 0
  PARAMETERS MATOADJUST=1
  PARAMETERS COMBINE=EQUI
  PARAMETERS ZONEMSG= 100
;-----------------------------------------------------------------------------;
  PHASE=LINKREAD
    LINKCLASS = LI.CAPCLASS             ; set linkclass variable to CAPCLASS
    LW.FFT=LI.DISTANCE/LI.SPEED*60      ; FREEFLOW TIME IN MINUTES
  ENDPHASE
;-----------------------------------------------------------------------------;
    PHASE=ADJUST
        FUNCTION v=vol[1]
        READ FILE=..\Output\?_TURNS.DAT
        ,T=TURN[1]

      ;  
      ; Conical function: T0 * ({2-b} - {a}*(1-x) + sqrt({a^2}*(1-x)^2 + {b^2}), where a=slope at capacity, b=(2a-1)/(2a-2)
        FUNCTION  TC[1] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Freeways
        FUNCTION  TC[2] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Rural Road (Expressway Curve)
        FUNCTION  TC[3] = T0 * (0.75  - 3  *(1 - (V/C)) + sqrt( 9  *(1 - (V/C))*(1 - (V/C))+1.5625)) ; Arterials
        FUNCTION  TC[4] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Minor Arterial
        FUNCTION  TC[5] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Collectors
        FUNCTION  TC[6] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[7] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[10] = T0                            ; Centroid Connectors (free-flow speeds)
    ENDPHASE
;-----------------------------------------------------------------------------;
  PHASE=ILOOP

    PATHLOAD PATH=TIME,PENI=1, VOL[1]=MI.1.TOTAL,
                               VOL[2]=MI.1.HBW,
                               VOL[3]=MI.1.HRS,
                               VOL[4]=MI.1.HSH,
                               VOL[5]=MI.1.HSR,
                               VOL[6]=MI.1.HSC,
                               VOL[7]=MI.1.NHB,
                               VOL[8]=MI.1.XX
  ENDPHASE
;-----------------------------------------------------------------------------;
ENDRUN

;==============================================================================
;==============================================================================
;  MIDDAY ASSIGNMENT
;==============================================================================
;==============================================================================
RUN PGM=HWYLOAD

NETI=..\Output\?base.net
FILEI TURNPENI = ..\Input\TURNPEN.DAT, MISSINGLINK=1
FILEI MATI[1] = ..\Output\MD_OD_ByPurp.MAT
FILEO NETO = ..\Output\?_MD.net

;-----------------------------------------------------------------------------;
  ; CAPACITY IN VPH/LANE:
  SPDCAP CAPACITY[1]=2000       ; FREEWAY
  SPDCAP CAPACITY[2]=1000       ; RURAL ROAD
  SPDCAP CAPACITY[3]=800        ; MAJOR ARTERIAL
  SPDCAP CAPACITY[4]=700        ; MINOR ARTERIAL
  SPDCAP CAPACITY[5]=600        ; COLLECTOR
  SPDCAP CAPACITY[6]=1500       ; LINEAR RAMP
  SPDCAP CAPACITY[7]=1200       ; LOOP RAMP
  SPDCAP CAPACITY[10]=9999      ; CENTROID CONNECTOR
;  SPDCAP CAPACITY[10]=3000      ; CENTROID CONNECTOR
;-----------------------------------------------------------------------------;
  PARAMETERS CAPFAC=4.84  ;MD 5-hour intensity factor
  PARAMETERS MAXITERS=70
  parameters relativegap=0.0002
  RMSE=0
  PARAMETERS GAP = 0
  PARAMETERS RAAD = 0
  PARAMETERS AAD = 0
  PARAMETERS MATOADJUST=1
  PARAMETERS COMBINE=EQUI
  PARAMETERS ZONEMSG= 100
;-----------------------------------------------------------------------------;
  PHASE=LINKREAD
    LINKCLASS = LI.CAPCLASS             ; set linkclass variable to CAPCLASS
    LW.FFT=LI.DISTANCE/LI.SPEED*60      ; FREEFLOW TIME IN MINUTES
  ENDPHASE
;-----------------------------------------------------------------------------;
    PHASE=ADJUST
        FUNCTION v=vol[1]
        READ FILE=..\Output\?_TURNS.DAT
        ,T=TURN[1]

      ;  
      ; Conical function: T0 * ({2-b} - {a}*(1-x) + sqrt({a^2}*(1-x)^2 + {b^2}), where a=slope at capacity, b=(2a-1)/(2a-2)
        FUNCTION  TC[1] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Freeways
        FUNCTION  TC[2] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Rural Road (Expressway Curve)
        FUNCTION  TC[3] = T0 * (0.75  - 3  *(1 - (V/C)) + sqrt( 9  *(1 - (V/C))*(1 - (V/C))+1.5625)) ; Arterials
        FUNCTION  TC[4] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Minor Arterial
        FUNCTION  TC[5] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Collectors
        FUNCTION  TC[6] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[7] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[10] = T0                            ; Centroid Connectors (free-flow speeds)
    ENDPHASE
;-----------------------------------------------------------------------------;
  PHASE=ILOOP

    PATHLOAD PATH=TIME,PENI=1, VOL[1]=MI.1.TOTAL,
                               VOL[2]=MI.1.HBW,
                               VOL[3]=MI.1.HRS,
                               VOL[4]=MI.1.HSH,
                               VOL[5]=MI.1.HSR,
                               VOL[6]=MI.1.HSC,
                               VOL[7]=MI.1.NHB,
                               VOL[8]=MI.1.XX

  ENDPHASE
;-----------------------------------------------------------------------------;
ENDRUN

;==============================================================================
;==============================================================================
;  PM PEAK PERIOD ASSIGNMENT
;==============================================================================
;==============================================================================
RUN PGM=HWYLOAD

NETI=..\Output\?base.net
FILEI TURNPENI = ..\Input\TURNPEN.DAT, MISSINGLINK=1
FILEI MATI[1] = ..\Output\P3_OD_ByPurp.MAT
FILEO NETO = ..\Output\?_P3.net

TURNVOLO[1]=..\Output\?_P3_turn.dbf, FORMAT=DBF, DEC=0

;-----------------------------------------------------------------------------;
  ; CAPACITY IN VPH/LANE:
  SPDCAP CAPACITY[1]=2000       ; FREEWAY
  SPDCAP CAPACITY[2]=1000       ; RURAL ROAD
  SPDCAP CAPACITY[3]=800        ; MAJOR ARTERIAL
  SPDCAP CAPACITY[4]=700        ; MINOR ARTERIAL
  SPDCAP CAPACITY[5]=600        ; COLLECTOR
  SPDCAP CAPACITY[6]=1500       ; LINEAR RAMP
  SPDCAP CAPACITY[7]=1200       ; LOOP RAMP
  SPDCAP CAPACITY[10]=9999      ; CENTROID CONNECTOR
;  SPDCAP CAPACITY[10]=3000      ; CENTROID CONNECTOR
;-----------------------------------------------------------------------------;
  PARAMETERS CAPFAC=2.96  ;3-hour intensity factor
  PARAMETERS MAXITERS=70
  parameters relativegap=0.0002
  RMSE=0
  PARAMETERS GAP = 0
  PARAMETERS RAAD = 0
  PARAMETERS AAD = 0
  PARAMETERS MATOADJUST=1
  PARAMETERS COMBINE=EQUI
  PARAMETERS ZONEMSG= 100
;-----------------------------------------------------------------------------;
  PHASE=LINKREAD
    LINKCLASS = LI.CAPCLASS             ; set linkclass variable to CAPCLASS
    LW.FFT=LI.DISTANCE/LI.SPEED*60      ; FREEFLOW TIME IN MINUTES
  ENDPHASE
;-----------------------------------------------------------------------------;
    PHASE=ADJUST
        FUNCTION v=vol[1]
        READ FILE=..\Output\?_TURNS.DAT
        ,T=TURN[1]

      ;  
      ; Conical function: T0 * ({2-b} - {a}*(1-x) + sqrt({a^2}*(1-x)^2 + {b^2}), where a=slope at capacity, b=(2a-1)/(2a-2)
        FUNCTION  TC[1] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Freeways
        FUNCTION  TC[2] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Rural Road (Expressway Curve)
        FUNCTION  TC[3] = T0 * (0.75  - 3  *(1 - (V/C)) + sqrt( 9  *(1 - (V/C))*(1 - (V/C))+1.5625)) ; Arterials
        FUNCTION  TC[4] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Minor Arterial
        FUNCTION  TC[5] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Collectors
        FUNCTION  TC[6] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[7] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[10] = T0                            ; Centroid Connectors (free-flow speeds)
    ENDPHASE
;-----------------------------------------------------------------------------;
  PHASE=ILOOP

    PATHLOAD PATH=TIME,PENI=1, VOL[1]=MI.1.TOTAL,
                               VOL[2]=MI.1.HBW,
                               VOL[3]=MI.1.HRS,
                               VOL[4]=MI.1.HSH,
                               VOL[5]=MI.1.HSR,
                               VOL[6]=MI.1.HSC,
                               VOL[7]=MI.1.NHB,
                               VOL[8]=MI.1.XX

  ENDPHASE
;-----------------------------------------------------------------------------;
ENDRUN

;==============================================================================
;==============================================================================
;  EVENING PERIOD ASSIGNMENT
;==============================================================================
;==============================================================================
RUN PGM=HWYLOAD

NETI=..\Output\?base.net
FILEI TURNPENI = ..\Input\TURNPEN.DAT, MISSINGLINK=1
FILEI MATI[1] = ..\Output\EV_OD_ByPurp.MAT
FILEO NETO = ..\Output\?_EV.net


;-----------------------------------------------------------------------------;
  ; CAPACITY IN VPH/LANE:
  SPDCAP CAPACITY[1]=2000       ; FREEWAY
  SPDCAP CAPACITY[2]=1000       ; RURAL ROAD
  SPDCAP CAPACITY[3]=800        ; MAJOR ARTERIAL
  SPDCAP CAPACITY[4]=700        ; MINOR ARTERIAL
  SPDCAP CAPACITY[5]=600        ; COLLECTOR
  SPDCAP CAPACITY[6]=1500       ; LINEAR RAMP
  SPDCAP CAPACITY[7]=1200       ; LOOP RAMP
  SPDCAP CAPACITY[10]=9999      ; CENTROID CONNECTOR
;  SPDCAP CAPACITY[10]=3000      ; CENTROID CONNECTOR
;-----------------------------------------------------------------------------;
  PARAMETERS CAPFAC=5.17  ;13-hour EV intensity factor
  PARAMETERS MAXITERS=70
  parameters relativegap=0.0002
  RMSE=0
  PARAMETERS GAP = 0
  PARAMETERS RAAD = 0
  PARAMETERS AAD = 0
  PARAMETERS MATOADJUST=1
  PARAMETERS COMBINE=EQUI
  PARAMETERS ZONEMSG= 100
;-----------------------------------------------------------------------------;
  PHASE=LINKREAD
    LINKCLASS = LI.CAPCLASS             ; set linkclass variable to CAPCLASS
    LW.FFT=LI.DISTANCE/LI.SPEED*60      ; FREEFLOW TIME IN MINUTES
  ENDPHASE
;-----------------------------------------------------------------------------;
    PHASE=ADJUST
        FUNCTION v=vol[1]
        READ FILE=..\Output\?_TURNS.DAT
        ,T=TURN[1]

      ;  
      ; Conical function: T0 * ({2-b} - {a}*(1-x) + sqrt({a^2}*(1-x)^2 + {b^2}), where a=slope at capacity, b=(2a-1)/(2a-2)
        FUNCTION  TC[1] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Freeways
        FUNCTION  TC[2] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Rural Road (Expressway Curve)
        FUNCTION  TC[3] = T0 * (0.75  - 3  *(1 - (V/C)) + sqrt( 9  *(1 - (V/C))*(1 - (V/C))+1.5625)) ; Arterials
        FUNCTION  TC[4] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Minor Arterial
        FUNCTION  TC[5] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Collectors
        FUNCTION  TC[6] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[7] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[10] = T0                            ; Centroid Connectors (free-flow speeds)
    ENDPHASE
;-----------------------------------------------------------------------------;
  PHASE=ILOOP

    PATHLOAD PATH=TIME,PENI=1, VOL[1]=MI.1.TOTAL,
                               VOL[2]=MI.1.HBW,
                               VOL[3]=MI.1.HRS,
                               VOL[4]=MI.1.HSH,
                               VOL[5]=MI.1.HSR,
                               VOL[6]=MI.1.HSC,
                               VOL[7]=MI.1.NHB,
                               VOL[8]=MI.1.XX

  ENDPHASE
;-----------------------------------------------------------------------------;
ENDRUN

;==============================================================================
;==============================================================================
;  AM PEAK HOUR ASSIGNMENT
;==============================================================================
;==============================================================================
RUN PGM=HWYLOAD

NETI=..\Output\?base.net
FILEI TURNPENI = ..\Input\TURNPEN.DAT, MISSINGLINK=1
FILEI MATI[1] = ..\Output\A1_OD_ByPurp.MAT
FILEO NETO = ..\Output\?_A1.net

TURNVOLO[1]=..\Output\?_A1_turn.dbf, FORMAT=DBF, DEC=0

;-----------------------------------------------------------------------------;
  ; CAPACITY IN VPH/LANE:
  SPDCAP CAPACITY[1]=2000       ; FREEWAY
  SPDCAP CAPACITY[2]=1000       ; RURAL ROAD
  SPDCAP CAPACITY[3]=800        ; MAJOR ARTERIAL
  SPDCAP CAPACITY[4]=700        ; MINOR ARTERIAL
  SPDCAP CAPACITY[5]=600        ; COLLECTOR
  SPDCAP CAPACITY[6]=1500       ; LINEAR RAMP
  SPDCAP CAPACITY[7]=1200       ; LOOP RAMP
  SPDCAP CAPACITY[10]=9999      ; CENTROID CONNECTOR
;  SPDCAP CAPACITY[10]=3000      ; CENTROID CONNECTOR
;-----------------------------------------------------------------------------;
  PARAMETERS CAPFAC=1.00  ;1-hour intensity factor
  PARAMETERS MAXITERS=70
  parameters relativegap=0.0002
  RMSE=0
  PARAMETERS GAP = 0
  PARAMETERS RAAD = 0
  PARAMETERS AAD = 0
  PARAMETERS MATOADJUST=1
  PARAMETERS COMBINE=EQUI
  PARAMETERS ZONEMSG= 100
;-----------------------------------------------------------------------------;
  PHASE=LINKREAD
    LINKCLASS = LI.CAPCLASS             ; set linkclass variable to CAPCLASS
    LW.FFT=LI.DISTANCE/LI.SPEED*60      ; FREEFLOW TIME IN MINUTES
  ENDPHASE
;-----------------------------------------------------------------------------;
    PHASE=ADJUST
        FUNCTION v=vol[1]
        READ FILE=..\Output\?_TURNS.DAT
        ,T=TURN[1]

      ;  
      ; Conical function: T0 * ({2-b} - {a}*(1-x) + sqrt({a^2}*(1-x)^2 + {b^2}), where a=slope at capacity, b=(2a-1)/(2a-2)
        FUNCTION  TC[1] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Freeways
        FUNCTION  TC[2] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Rural Road (Expressway Curve)
        FUNCTION  TC[3] = T0 * (0.75  - 3  *(1 - (V/C)) + sqrt( 9  *(1 - (V/C))*(1 - (V/C))+1.5625)) ; Arterials
        FUNCTION  TC[4] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Minor Arterial
        FUNCTION  TC[5] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Collectors
        FUNCTION  TC[6] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[7] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[10] = T0                            ; Centroid Connectors (free-flow speeds)
    ENDPHASE
;-----------------------------------------------------------------------------;
  PHASE=ILOOP

    PATHLOAD PATH=TIME,PENI=1, VOL[1]=MI.1.TOTAL,
                               VOL[2]=MI.1.HBW,
                               VOL[3]=MI.1.HRS,
                               VOL[4]=MI.1.HSH,
                               VOL[5]=MI.1.HSR,
                               VOL[6]=MI.1.HSC,
                               VOL[7]=MI.1.NHB,
                               VOL[8]=MI.1.XX

  ENDPHASE
;-----------------------------------------------------------------------------;
ENDRUN

;==============================================================================
;==============================================================================
;  PM PEAK HOUR ASSIGNMENT
;==============================================================================
;==============================================================================
RUN PGM=HWYLOAD

NETI=..\Output\?base.net
FILEI TURNPENI = ..\Input\TURNPEN.DAT, MISSINGLINK=1
FILEI MATI[1] = ..\Output\P1_OD_ByPurp.MAT
FILEO NETO = ..\Output\?_P1.net

TURNVOLO[1]=..\Output\?_P1_turn.dbf, FORMAT=DBF, DEC=0

;-----------------------------------------------------------------------------;
  ; CAPACITY IN VPH/LANE:
  SPDCAP CAPACITY[1]=2000       ; FREEWAY
  SPDCAP CAPACITY[2]=1000       ; RURAL ROAD
  SPDCAP CAPACITY[3]=800        ; MAJOR ARTERIAL
  SPDCAP CAPACITY[4]=700        ; MINOR ARTERIAL
  SPDCAP CAPACITY[5]=600        ; COLLECTOR
  SPDCAP CAPACITY[6]=1500       ; LINEAR RAMP
  SPDCAP CAPACITY[7]=1200       ; LOOP RAMP
  SPDCAP CAPACITY[10]=9999      ; CENTROID CONNECTOR
;  SPDCAP CAPACITY[10]=3000      ; CENTROID CONNECTOR
;-----------------------------------------------------------------------------;
  PARAMETERS CAPFAC=1.00  ;1-hour intensity factor
  PARAMETERS MAXITERS=70
  parameters relativegap=0.0002
  RMSE=0
  PARAMETERS GAP = 0
  PARAMETERS RAAD = 0
  PARAMETERS AAD = 0
  PARAMETERS MATOADJUST=1
  PARAMETERS COMBINE=EQUI
  PARAMETERS ZONEMSG= 100
;-----------------------------------------------------------------------------;
  PHASE=LINKREAD
    LINKCLASS = LI.CAPCLASS             ; set linkclass variable to CAPCLASS
    LW.FFT=LI.DISTANCE/LI.SPEED*60      ; FREEFLOW TIME IN MINUTES
  ENDPHASE
;-----------------------------------------------------------------------------;
    PHASE=ADJUST
        FUNCTION v=vol[1]
        READ FILE=..\Output\?_TURNS.DAT
        ,T=TURN[1]

      ;  
      ; Conical function: T0 * ({2-b} - {a}*(1-x) + sqrt({a^2}*(1-x)^2 + {b^2}), where a=slope at capacity, b=(2a-1)/(2a-2)
        FUNCTION  TC[1] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Freeways
        FUNCTION  TC[2] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Rural Road (Expressway Curve)
        FUNCTION  TC[3] = T0 * (0.75  - 3  *(1 - (V/C)) + sqrt( 9  *(1 - (V/C))*(1 - (V/C))+1.5625)) ; Arterials
        FUNCTION  TC[4] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Minor Arterial
        FUNCTION  TC[5] = T0 * (0.667 - 2.5*(1 - (V/C)) + sqrt(6.25*(1 - (V/C))*(1 - (V/C))+1.7778)) ; Collectors
        FUNCTION  TC[6] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[7] = T0 * (0.875 - 5  *(1 - (V/C)) + sqrt(25  *(1 - (V/C))*(1 - (V/C))+1.2657)) ; Slip and Loop Ramps (Expressway Curve)
        FUNCTION  TC[10] = T0                            ; Centroid Connectors (free-flow speeds)
    ENDPHASE
;-----------------------------------------------------------------------------;
  PHASE=ILOOP

    PATHLOAD PATH=TIME,PENI=1, VOL[1]=MI.1.TOTAL,
                               VOL[2]=MI.1.HBW,
                               VOL[3]=MI.1.HRS,
                               VOL[4]=MI.1.HSH,
                               VOL[5]=MI.1.HSR,
                               VOL[6]=MI.1.HSC,
                               VOL[7]=MI.1.NHB,
                               VOL[8]=MI.1.XX

  ENDPHASE
;-----------------------------------------------------------------------------;
ENDRUN

;==============================================================================
;==============================================================================
;  DAILY TOTAL
;==============================================================================
;==============================================================================
run pgm=network
; Combine period volumes for daily total
filei
   neti[1]= ..\Output\?base.net
   neti[2]= ..\Output\?_a3.net
   neti[3]= ..\Output\?_md.net
   neti[4]= ..\Output\?_p3.net
   neti[5]= ..\Output\?_ev.net
   neti[6]= ..\Output\?_a1.net
   neti[7]= ..\Output\?_p1.net
   

fileo
   neto= ..\Output\?_day.net

merge record=false 

   A3VOL = li.2.V_1
   P3VOL = li.4.V_1
   MDVOL = li.3.V_1
   EVVOL = li.5.V_1
   A1VOL = li.6.V_1
   P1VOL = li.7.V_1
   DAYVOL= li.2.V_1 + li.3.V_1 + li.4.V_1 +li.5.V_1
   DAYVOLT= li.2.VT_1 + li.3.VT_1 + li.4.VT_1 + li.5.VT_1
   maxvc = max(li.2.vc_1, li.3.vc_1, li.4.vc_1, li.5.vc_1)
   dayvmt = li.2.vdt_1 + li.3.vdt_1 + li.4.vdt_1 + li.5.vdt_1
   dayvht = li.2.vht_1 + li.3.vht_1 + li.4.vht_1 + li.5.vht_1
   if (dayvht > 0)
       avgspd = dayvmt / dayvht
   else
       avgspd = 0
   endif
 if (_outcount=0)   
   print file=..\Output\?_day.csv, csv=T, form=(6.0), 
      list='A,B,capclass,lanes,speed,cnt_A1,cnt_P1,cnt_A2,cnt_P2,cnt_A3,cnt_P3,cnt_MD,cnt_EV,cnt_Day,cnt_DayT,A1vol,P1vol,A3vol,P3vol,MDvol,EVvol,DAYVOL,DAYVOLT,index,MaxVC'
   _outcount=1
 endif

   print file=..\Output\?_day.csv, csv=T, form=(6.0), 
      list=A,B,capclass,lanes,speed,cnt_A1,cnt_P1,cnt_A2,cnt_P2,cnt_A3,cnt_P3,cnt_MD,cnt_EV,cnt_Day,cnt_DayT,A1vol,P1vol,A3vol,P3vol,MDvol,EVvol,DAYVOL,DAYVOLT,index,MaxVC(4.2)
   
endrun