ACES PSC Module:  Equations File   (Last modified:  10/7/09)

 This file contains a list of all equations used by the analysis routine.
 When performing an analysis calculations are performed in the strict
 order specified within this file.   New parameters can be defined and
 incorporated into the this Equations file but they must first be defined
 in the Parameters List file.   

 EDITING INSTRUCTIONS

 The rules for editing this file and creating and defining equations are
 given below.

 GENERAL RULES

 Since this module mimics a "spreadsheet" equations must be inserted
 into the list in the strict order in which you require them to be
 executed when performing a design.  Consequently they must conform
 to the following rules:

  (1)  Variables must first be defined in the Parameters List file.
  (2)  Blank lines may be used anywhere to delineate equation blocks.
  (3)  Comment lines must be prefaced with the exclamation symbol (!)
  (4)  Each equation may be up to 120 characters long.
  (5)  Allowable operands are ‘+’, ‘-‘, ‘*’, ‘/’ and ‘^’ (denotes exponential)
  (6)  Variables must appear in equations without curly braces.
  (7)  Undefined variables will be allocated the value ‘1’.
  (8)  Any number of bracket groups may be used in an equation.
  (9)  Nested brackets are NOT permitted.
 (10) Spaces in the equation will be ignored.
 (11) Constants may be specified as decimal, integer or exponential
 (12) String (text) variables must contain the '$' symbol (e.g. Note$)
 (13) String variables can only be assigned a maximum of 44 characters
 (14) Comments on an IF statement line must not include the ( ) symbols
 (15) Only equations that lie within the **** boundaries will be used

 ________________________________________________________________________

 DISPLAYING EQUATIONS in the CALCULATION LOG

    .  Equations have been listed sequentially in the equation list to
       follow the logical order in which the design is usually performed.
       Blocks of equations therefore correspond to each tab in the 
       module.  The tabs have been internally numbered 1-9 and the 
       first block of equations in any tab is also given that same number
       e.g.               !4   PRESTRESS LOSSES TAB

    .  When the <Check calcs> button at the bottom of a tab is clicked 
       ACES will display the calculation log of all equations associated
       with that tab.

    .  Sub-groups of equations represent various frames and calculation
       buttons within a given tab.  They have been given a label and, 
       where applicable, a letter designating that sub-group  e.g.,
                !7b  SHEAR DESIGN - Flexure shear cracking

    .  Only the equations relating to this block will be displayed in the
       calculation log when the 'View calcs' button on that dialog box (or
       in that frame) is clicked.

    .  When the <Full calcs> button at the bottom of a tab is clicked
       ACES will display the full calculation log i.e., all equations
       associated with the entire analysis.

    .  A block or sub-block is ended when another tab group number or
       sub-group number is encountered.

    .  !0 signifies that the equations in that block are not to be associated
       with any button.

    .  !+ signifies that all equations in that block will always be displayed
       in the log.

    .  !- signifies that equations in that block will not be displayed in the
       log for subgroups.

 ________________________________________________________________________

 GENERAL FUNCTIONS

    .  The MAX and MIN functions 

       The MAXimising and MINimising functions operate on variables 
       and numbers only - no expressions are allowed.  Variables and
       numbers may be signed  (e.g.  -ve) and up to 20 values only are
       allowed within the brackets.  Note that the ABS and TAN 
       functions CANNOT be included in a MAX/MIN list.

    .  TAN (tangent) function.  This takes the form:
                              var = TAN (Angle)
       'Angle' must be specified in radians.  If it is not already in radians
       then convert it using one of the nine temporary work variables
       (WorkVarX) described in the Parameters List file and the two
       expressions:
                              WorkVarX = (Angle*3.14159/180.)
                              var = TAN (WorkVarX)
       Note that the ABS and TAN functions CANNOT be included together
       on the same line nor in a MAX/MIN list.  In such cases a
       temporary parameter (WorkVarX) must be used to calculate an
       intermediate value.

    .  ABS (absolute value) function.  This takes the form:
                              var = ABS (xxx)
       where 'var' is any parameter in the parameter list and 'xxx' is either
       an optionally signed variable or a number.  It MUST be included
       within normal brackets (as shown).   Note that the ABS and TAN
       functions CANNOT be included together on the same line nor in a
       MAX/MIN list.  In such cases a temporary parameter (WorkVarX)
       must be used to calculate an intermediate value.  Note also that in
       a statement such as var = -ABS (xxx)  var will always be negative.

    .  OUTPUT function

       The function is specified as:  OUTPUT a, b, c, d, .... etc  where a,
       b, c, etc  are parameters that have already been defined and used
       in the equation list.  ACES will simply display the specified values,
       separated by commas.   This can be very useful in LOOP and
       ITERATE functions to track the behaviour of one or more variables.
       A maximum of 20 parameters can be listed in OUTPUT ......

       Note that the display of values will occur even if the NOTRACK 
       flag has been invoked (refer to the LOOP function below for a
       description of the NOTRACK flag).


    .  IF function.  There are two types of IF function:

       (1) Single line form:  The function can be expressed either as:
                       IF (xxx > yyy)  expression
                                       or as:
                       IF (xxx > yyy) THEN expression

       (2) Multiple line form.  This must be expressed as:
                       IF (xxx > yyy) THEN
                            expressions1
                       ELSE
                            expressions2
                       END IF

            Both expressions1 and expressions2 may contain up to 20
            statements (or none at all) and the operands can only be:
            ">", "<" or "=".   Values 'xxx' and 'yyy' can only be variables or
            numbers and must be included within normal brackets.
            Comments on an IF statement line must not include embedded
            quotation ( ) symbols.


    .  LOOP function.  This takes the form:

            LOOP  i = a, b NOTRACK
                  expression 1
                  function 2
	   ...... etc ......
            END LOOP

       NOTRACK is an optional attribute - it suppresses the tracking and
       display of intermediate results during the looping process.

       If any parameter name within 'expression' contains the character 'i'
       then the 'i' itself will be replaced with the number corresponding to
       the loop variable 'i'.

       For example, in the calculation of tension forces in PS strands the
       general equation for setting the bar depth (Dpbari) of bars in row i
       is:  “Dpbari = D – Ybarri”.  Using the LOOP function, this can be
       performed for bar rows 1 to 8 using three lines of code  viz:

             LOOP  i = 1,8
                   Dpbarri = D - Ybarri
             END LOOP

       Internally ACES will create 8 equations of the form:
       Dpbarr1 = D - Ybarr1, Dpbarr2 = D - Ybarr2, etc.

       A maximum of 20 expressions and/or functions can be included in
       the loop and the loop variables 'a' and 'b' can only lie betweeb 0 and
       16  (e.g., LOOP i = 1, 20 would not be accepted by the system).


    .  ITERATE funtion.  This must take the form:

        ITERATE
          .........
        END ITERATE

              or

        ITERATE ON param IN STEPS OF xx  NOTRACK
            ......  (up to twenty statements allowed) .....
        END ITERATION WHEN (param ><= number) OR AFTER n CYCLES

        where:  'param' = any parameter in the parameter list

                         'xx' = the number to be added to 'param' at the
                                 beginning of each cycle. This may be zero if
                                 the iteration is naturally convergent.

                         ' n' = maximum number of cycles to perform if the
                                 iteration is not converging quickly enough.

       The quantity evaluated, (param ><= number) must be enclosed in
       normal brackets.  Note that the tolerance for equality has been
       set to 1%.  That is, if two numbers are within 1% of each other,
       they are considered equal. 'number' may be a parameter or a 
       number.

       NOTRACK is optional and suppresses the tracking and display of
       intermediate results during the iteration process.

       Example:

         ITERATE ON WorkVar1 IN STEPS OF 5
                WorkVar2 = WorkVar1 * 2.
         END ITERATION WHEN (WorkVar1 = 25) OR AFTER 30 CYCLES

 ________________________________________________________________________

SPECIAL FUNCTIONS

    .  SORT functions:
                                 SORT_LOW_HIGH (p1i,p2i,r1i,r2i) NOTRACK
                                 SORT_HIGH_LOW (p1i,p2i,r1i,r2i) NOTRACK

       These two functions sort a set of 8 values 'p1i' in ascending and/or
       descending order respectively and place the sorted results into the
       set labelled 'r1i'.  Values in set 'p2i' correspond to those in set 'p1i'
       and are placed into 'r2i' during the sort process.  In effect the
       function sorts the set 'p1i' and its corresponding values 'p2i' and
       stores the results into 'r1i' and 'r2i'.  Note, however, that 'p1i' and
       'p2i' must contain 8 values each.

       NOTRACK suppresses the tracking and display of sorted output.
       The functions have been set up specifically for sorting strand and
       passive reinforcement bars and their corresponding forces in order
       to determine if they lie in the tension zone of the section.

       Example:      SORT_LOW_HIGH (Ybarri, Fbari, Dsorti, Fsorti)


    .  COMBINE and SORT functions:

              COMBINE_SORT_LOW_HIGH (p1i, p2i, p3i, p4i, r1i, r2i, r3i)
              COMBINE_SORT_HIGH_LOW (p1i, p2i, p3i, p4i, r1i, r2i, r3i)

       This is similar to the SORT function described above with the
       exception that the 'p1i' and 'p3i' data sets and their corresponding
      'p2i' and 'p4i' sets are first combined before sorting.  The 8 values in
       each of the individual 'p1i' and 'p3i' sets then appear as 16 sorted
       values in the final 'r1i' set.  This function is used to combine the
       strand table ('p1i'/'p3i') and the passive reinforcement tables into
       one consolidated table which is then sorted by bar distance ('r1i')
       and the corresponding bar forces ('r2i').

       The values in 'r3i' are flags that indicate which set the 'r1i' values
       originally came from (either 'p1i' or 'p3i').  A value of 1 indicates 'p1i'
       and 2 indicates 'p3i'.  This information is required by the strain
       compatibility routine to differentiate a PS strand from a reinforcing 
       bar when calculating the yield strain.

       Example:

       COMBINE_SORT_HIGH_LOW (Dbari,Fbari,Dpfi,Fbfi,Dsti,Fsti,Typei)


    .  LOCATE Ultimate Neutral Axis Depth (Dcb) Functions:

                     LOCATE_DCB_USING_SIMPLIFIED
                                           and
                     LOCATE_DCB_USING_STRAIN_COMPATIBILITY

       Both functions use the values returned by the COMBINE_SORT
       functions.  They locate the neutral axis depth (Dcb) and populate
       the bar force table automatically.  Although these functions could
       have been replaced by inserting relevant equations directly into the
       equations file it has been found that the hard-coded routines
       implemented here-in execute much faster.

       The '_SIMPLIFIED' method assumes:  (1) a concrete stress block of
       0.85f'c; and (2) that all bars and strands have yielded (except those
       in the compression zone, which are excluded from the analysis).

       The '_STRAIN _COMPATIBILITY' method derives the concrete 
       ultimate strain from the parameter 'ucu' which is editable on the
       'Ult M Check' tab but which defaults to 0.003.  The function 
       incrementally adjusts the depth Dcb until the concrete force 
       balances the steel forces.  Bars and strands in the compression 
       zone are excluded and the bar/strand forces are reduced if the
       bar/strand strain is below yield.  (This is commonly refered to as
       'Partial Prestressing').


***************************** Start of Equations *****************************

!-  CALCULATION OF DEFAULT VALUES

IF (f`cg = 0.) THEN f`cg = 50.
IF (Eg = 0.) THEN Eg = 0.043 * (2500.^1.5) * (f`cg^0.5)
IF (f`cs = 0.) THEN f`cs = 40.
IF (Es = 0.) THEN Es = 0.043 * (2400.^1.5) * (f`cs^0.5)
IF (f`cmi = 0.) THEN f`cmi = 40.
IF (Egmi = 0.) THEN Egmi = 0.043 * (2500.^1.5) * (f`cmi^0.5)
IF (f`cmt = 0.) THEN f`cmt = 30.
IF (Egmt = 0.) THEN Egmt = 0.043 * (2500.^1.5) * (f`cmt^0.5)
IF (ø = 0.) THEN ø = 0.9
IF (O = 0.) THEN O = 0.9
IF (øm= 0.) THEN øm= 0.9
IF (Om= 0.) THEN Om= 0.9
IF (ucu = 0.) THEN ucu = 0.003

!2   SECTION PROPERTIES TAB - Section property calculations

Zt = Ig / (Dg - Yb)
Zb = Ig / Yb

IF (Ws = 0.) THEN Ws = Wtf
IF (Es = 0.) THEN Es = Eg

n =  Es/Eg

IF (CalcProp = 1) THEN
   CalcPrp$ = Composite section props recalculated by ACES
   D = Dg + Ts
   As = Ws * Ts
   Ys = Dg + Ts/2.
   Ac = Ag + As*n
   Yc = (Ag * Yb + As*n*Ys) / (Ag + As*n)
   WorkVar1 = Yc - Yb
   WorkVar2 = Ys - Yc
   Ic = Ig + Ag*WorkVar1^2 + n*Ws*Ts^3/12. + n*As*WorkVar2^2
   WorkVar1 = D - Yc
   Zst = Ic / (WorkVar1*n)
   WorkVar1 = Dg - Yc
   Zsb = Ic / (WorkVar1*n)
   Zgt = Ic / (WorkVar1)
   Zgb = Ic / Yc
ELSE
   CalcPrp$ = Composite props NOT recalculated by ACES 
END IF

!  Calculation of shear flow at centroid and web/flange

IF (CalcProp = 1) THEN
   Ybflan = (Yc - Ybf)                     ! Bottom flange to composite centroid
   Hwebfl = (Yc - Dwf)                     ! Web/flange to composite centroid
   Ywebfl = Hwebfl / 2.                    ! Centroidal offset of web block
   Areac = Abf + Hwebfl*Bw                 ! Area of bottom flange + web block
   Ymomsc = Abf*Ybflan + Hwebfl*Bw*Ywebfl  ! Moment of area of section below centroid
   Ybarc = Ymomsc / Areac                  ! Centroidal offset of section below centroid
   Qna = Areac * Ybarc                     ! Shear flow at centroid
END IF

IF (Qna < 0.) THEN Qna = -Qna
IF (CalcProp = 1) THEN Qfw = Abf * Ybflan  ! Shear flow at web/flange interface
IF (Qfw < 0.) THEN Qfw = -Qfw

!2a   SECTION PROPERTIES TAB - Non-prestressed reinforcement

Nrfbarst = 0
  Arft     = 0
Ynbarsrf = 0

LOOP  i = 1,20
  Nrfbarst = Nrfbarst + Nrfbrsi
END LOOP

LOOP  i = 1,20
  WorkVar1 = 3.14159 * BarDtri * BarDtri / 4
  Arfi = WorkVar1 * Nrfbrsi
END LOOP

LOOP  i = 1,20
  Arft = Arft + Arfi
END LOOP

Arl = Arft          ! Total area of passive reinforcement

LOOP  i = 1,20
  YbxNrfi = Yrfi * Arfi
END LOOP

LOOP i = 1,20
  Ynbarsrf = Ynbarsrf + YbxNrfi
END LOOP


!3a PS ESTIMATES TAB - Differential Shrinkage

u = u2 - u1
n = Es/Eg
ec = Ys-Yc
Fshr = u*Es*As*1E-9*(1-2.71828^-Rcf)/Rcf   ! Calculation of shrinkage force
Mshr = Fshr*ec/1000
SumNA = n*As + Ag
Fsus = Fshr*1000*(1/As-n/SumNA)
Fsug = -Fshr*1000/SumNA

fts = Fsus - Mshr*1E6/Zst
fbs = Fsus - Mshr*1E6/Zsb
ftg = Fsug - Mshr*1E6/Zgt
fbg = Fsug + Mshr*1E6/Zgb

!3b PS ESTIMATES TAB - Prestress Forces

IF (DebondBs = 0) THEN
   Nbarb1 = Nbart1 - Nbard1
   Nbarb2 = Nbart2 - Nbard2
   Nbarb3 = Nbart3 - Nbard3
   Nbarb4 = Nbart4 - Nbard4
   Nbarb5 = Nbart5 - Nbard5
   Nbarb6 = Nbart6 - Nbard6
   Nbarb7 = Nbart7 - Nbard7
   Nbarb8 = Nbart8 - Nbard8
   IncDBar$ = Debonded bars are NOT included in analysis
ELSE
   Nbarb1 = Nbart1
   Nbarb2 = Nbart2
   Nbarb3 = Nbart3
   Nbarb4 = Nbart4
   Nbarb5 = Nbart5
   Nbarb6 = Nbart6
   Nbarb7 = Nbart7
   Nbarb8 = Nbart8
   IncDBar$ = Debonded bars are included in the analysis
END IF

YbxNb1 = Ybarr1 * Nbarb1
YbxNb2 = Ybarr2 * Nbarb2
YbxNb3 = Ybarr3 * Nbarb3
YbxNb4 = Ybarr4 * Nbarb4
YbxNb5 = Ybarr5 * Nbarb5
YbxNb6 = Ybarr6 * Nbarb6
YbxNb7 = Ybarr7 * Nbarb7
YbxNb8 = Ybarr8 * Nbarb8

!   Find total number of bars in layers 1 to 8 (Nbart1 to Nbart8)

Nbars = 0

LOOP  i = 1,8
  Nbars = Nbars + Nbarti
END LOOP

!   Find total number of debonded bars in layers 1 to 8

Ndbars = 0

LOOP  i = 1,8
  Ndbars= Ndbars+ Nbardi
END LOOP

!   Find total number of bars used in the analysis (may/may not include debonded bars)

Nbbars = 0

LOOP  i = 1,8
  Nbbars= Nbbars+ Nbarbi
END LOOP

!   Find total moment of area of bars in each layer

Ynbars = 0

LOOP  i = 1,8
  Ynbars= Ynbars+ YbxNbi
END LOOP

!   Calc Pj, centroid of bar group, eccentricity of group

Pj = Nbbars*Pult*Jf
Ycgs = Ynbars / Nbbars
e = Yb - Ycgs


!4  PRESTRESS LOSSES TAB - Steam Relaxation Losses

k5a = 1 + (Jf-0.7)*0.5/0.1
k5b = (Jf-0.4)/0.3
WorkVar3 = 0.
k5 = MAX (k5a,k5b,WorkVar3)

Lsrl = 0.1*k5/1.5
Lsr = 100*Lsrl
Prl = -Lsrl*Pj
Pjr = Pj + Prl

!   PRESTRESS LOSSES TAB - Elastic Deformation Losses

Ap = Nbbars * Aps
fcgs = -Pjr*1000*(1/Ag+e*e/Ig) + Msw*1E6*e/Ig
Pelastic = fcgs*Ep*Ap/(Egmt*1000)
Ledl = -Pelastic*100/Pj
Pt = Pjr + Pelastic
Ltr = Pt*100/Pj

!   PRESTRESS LOSSES TAB - Shrinkage Losses

Ac= n*As + Ag
Pshr = -us*Ep*Ap*1E-9/(1+15*Arl/Ac)
Lshr = -Pshr*100/Pj
Prs = Pt + Pshr

!   PRESTRESS LOSSES TAB - Creep loss parameters

At = Ag + Ws * Ts
th = 2*At/(Gp+0.5*Vp)
Fratio = f`cmt/f`cg

!4e PRESTRESS LOSSES TAB - Creep losses due to PS + SW

øcc = øccb*k2*k3
Occ = øcc              ! HTML Reports cannot handle ø symbol
Occb = øccb          ! HTML Reports cannot handle ø symbol
O = ø                    ! HTML Reports cannot handle ø symbol
fcscgs = -Pt*1000/Ag - Pt*1000*e*e/Ig + Msw*1E6*e/Ig
ucc1 = 1.E6*fcscgs*øcc/(Eg)

!4f PRESTRESS LOSSES TAB - Creep losses due to Deck + SDL

Fdeck = Mslab*1E6*(Yb-Ycgs)/Ig
Fsdl  = Msdl * 1E6*(Yc-Ycgs)/Ic
øcc2  = øccb * k2s*k3s
Occ2  = øcc2                                        ! HTML Reports cannot handle ø symbol
ucc2  = øcc2*1E6*(Fdeck+Fsdl)/Eg

!   PRESTRESS LOSSES TAB - Total creep losses

ucc    =  ucc1 + ucc2
Pcreep = ucc*Ep*Ap/1E9
Lcr    =  -Pcreep*100/Pj
P      =  Pt+Pshr+Pcreep
Ltt    =  P*100/Pj


!5  SERVICEABILITY CHECK - Girder stresses at transfer

f`csat = 0.5*f`cmt^0.5
f`csac = 0.6*f`cmt
f`at   = 0.5*f`cg^0.5
f`ac   = 0.6*f`cg
Mpte   = -Pt*e/1000
fgt1   = -Pt*1000/Ag
fgb1   =  fgt1
fgt2   = -Mpte*1E6/Zt
fgb2   =  Mpte*1E6/Zb
fgt3   = -Msw *1E6/Zt
fgb3   =  Msw*1E6/Zb
fgt4   =  fgt1+fgt2+fgt3
fgb4   =  fgb1+fgb2+fgb3

!5b SERVICEABILITY CHECK - Design moment stresses

Mpe    = -P*e/1000
fgtss1 = -P*1000/Ag
fgbss1 =  fgtss1
fgtss2 = -Mpe*1E6/Zt
fgbss2 =  Mpe*1E6/Zb
fgtss3 = -Msw*1E6/Zt
fgbss3 =  Msw*1E6/Zb
fgtss4 = -Mslab*1E6/Zt
fgbss4 =  Mslab*1E6/Zb
fstss5 = -Msdl*1E6*SDLf/Zst
fsbss5 = -Msdl*1E6*SDLf/Zsb
fgtss5 = -Msdl*1E6*SDLf/Zgt
fgbss5 =  Msdl*1E6*SDLf/Zgb
fstss6 =  fts
fsbss6 =  fbs
fgtss6 =  ftg
fgbss6 =  fbg
fstss7 = -Mll*1E6/Zst
fsbss7 = -Mll*1E6/Zsb
fgtss7 = -Mll*1E6/Zgt
fgbss7 =  Mll*1E6/Zgb
fstss8 = -Mhvl*1E6/Zst
fsbss8 = -Mhvl*1E6/Zsb
fgtss8 = -Mhvl*1E6/Zgt
fgbss8 =  Mhvl*1E6/Zgb

!   SERVICEABILITY CHECK - Final design DL+SDL+LL stresses

fstll  =  fstss5+fstss6+fstss7
fsbll  =  fsbss5+fsbss6+fsbss7
fgtll  =  fgtss1+fgtss2+fgtss3+fgtss4+fgtss5+fgtss6+fgtss7
fgbll  =  fgbss1+fgbss2+fgbss3+fgbss4+fgbss5+fgbss6+fgbss7

!   SERVICEABILITY CHECK - Final design DL+SDL+HLV stresses

fstsv  = fstss5+fstss6+fstss8
fsbsv  = fsbss5+fsbss6+fsbss8
fgtsv  = fgtss1+fgtss2+fgtss3+fgtss4+fgtss5+fgtss6+fgtss8
fgbsv  = fgbss1+fgbss2+fgbss3+fgbss4+fgbss5+fgbss6+fgbss8

!   SERVICEABILITY CHECK - MAXIMA

Msdlf = Msdl*SDLf 
fmaxsllt = fgbll

IF (fmaxsllt < 0.) THEN fmaxsllt = 0.

fmaxsllc = MIN (fstll,fsbll,fgtll)
fmaxssvt = fgbsv

IF (fmaxssvt < 0.) THEN fmaxssvt = 0.
fmaxssvc = MIN (fstsv,fsbsv,fgtsv)


!6b ULTIMATE MOMENT CHECK - Input data & constants 

Dp = D - Ycgs
fp = 1000. * Pult / Aps
WorkVar1 = Ap * fp
WorkVar2 = (Ast-Asc) * fsy
k2u = (WorkVar1 + WorkVar2) / (Ws*Dp*f`cs)
v = Srf - 0.007 * (f`cs-28)
fpu = fp * (1-k1u*k2u/v)
fuc = Srf * f`cs

!   ULTIMATE MOMENT CHECK - Tension forces

LOOP i = 1,8                       ! Strand depths
   Dpbari = D - Ybarri
END LOOP

LOOP i = 1,20                       ! Passive reinforcement depths
   Dprfi = D - Yrfi
END LOOP

LOOP i = 1,8                       ! Strand areas
   Apbari = Nbarbi * Aps
END LOOP

!   Passive reinforcement areas (Arfi) are already calculated

LOOP i = 1,8                         ! Forces in strand
   Fbari = Apbari * fpu/1000.
   IF (Fbari = 0.) THEN Dpbari = 0.
END LOOP

LOOP i = 1,20                         ! Forces in passive reo
   Fbrfi = Arfi * fsy / 1000.
   IF (Fbrfi = 0.) THEN Dprfi = 0.
END LOOP

!   Combine and sort prestressed and non-prestressed R/F
!   in descending order of distance from top fibre

 COMBINE_SORT_HIGH_LOW (Dpbari,Fbari,Dprfi,Fbrfi,Dsorti,Fsorti,Typei) NOTRACK

!   Calculate depth to neutral axis (Dcb) & steel/concrete capacities.
!   Use values stored in Dsorti and Fsorti and base capacity
!   calculations on the type of method selected (see below).
!
!  SIMPLIFIED METHOD - Assumptions
!
!  1. The ultimate concrete stress state is assumed to be
!     rectangular down to the neutral axis (depth = Dcb)
!
!  2. The ultimate concrete compressive stress Fuc = 0.85F’c
!
!  3. All steel (including passive reinforcement) below the neutral
!     axis is at yield.  Steel in compression is ignored.
!
!  4. The ultimate condition occurs when the compressive force in
!     the concrete balances the tensile forces in the steel.  The
!     moment capacity is then given by:
!              Mc = - Fuc*Ws*Dcb*Dcb/(2.*1000)
!
!
!  STRAIN-COMPATIBILITY METHOD - Assumptions
!
!  1. The ultimate concrete stress state is assumed to be
!     triangular down to the neutral axis (depth = Dcb).
!
!  2. The ultimate concrete strain at the top surface is given by
!     'ucu' (usually 0.003 microstrain).  The strain at each steel level
!     is determined by linear interpolation of a straight line running
!     from 'ucu' and through the neutral axis.
!
!  3. If the interpolated strain exceeds the yield strain for its
!     particular grade of steel, the full yield force for that bar is included.
!
!  4. If the interpolated strain is less that the yield strain, the
!     elastic force at that strain is used.
!
!  5. The ultimate condition occurs when the compressive force in
!     the concrete balances the tensile forces in the steel.  The
!     concrete stress block is triangular and the ultimate moment
!     capacity is given by:  Mc = - Fuc*Ws*Dcb*0.5*Dcb/(3.*1000)

IF (CalcDcb = 0) THEN
    McNote$ = Based on simplified method
    CalcDcb$ =Dcb calculated using simplified method
    LOCATE_DCB_USING_SIMPLIFIED_METHOD
ELSE
    McNote$ = Based on strain compatibility method
    CalcDcb$ =Dcb calc'd using strain compatibility method
    LOCATE_DCB_USING_STRAIN_COMPATIBILITY
END IF

!   Calculate area of passive reinforcement in tension zone

Ast = 0.

LOOP i = 1,20
   WorkVar1 = Arfi
   IF (Dprfi > Dcb) Ast = Ast + WorkVar1
END LOOP

LOOP i = 1,8
   IF (Fbari = 0.)  Yelds$i =
END LOOP

LOOP i = 1,20
   IF (Fbrfi = 0.)  Yeldr$i =
END LOOP

!   Calculate area of non-P/S steel in compression

Asc = Arft - Ast
Dcbnok = Dcb

!6c ULTIMATE MOMENT CHECK - Final Ultimate Design Moments

Musw = Msw * LFsw
Muslab = Mslab * LFslab
Musdl = Msdl * LFsdl
Mull = Mll * LFll
Muhvl = Mhvl * LFhvl
Mu1 = Musw + Muslab + Musdl + Mull
Mu2 = Musw + Muslab + Musdl + Muhvl
Mu = Mcs + Mc
øMu =  øm * Mu
Om = øm                            ! HTML report can`t handle ø symbol
OMu  = øMu                       ! HTML report can`t handle ø symbol
Mmax = MAX (Mu1,Mu2)

!7a  SHEAR DESIGN - Ultimate Design Parameters

WorkVar1 = LFsw*Vsw + LFslab*Vslab + LFsdl*Vsdl + LFll*Vll

IF (Vult = 0) THEN Vult = WorkVar1

st = 0.33*f`cg^0.5
f`cf = 0.6*f`cg^0.5
r = (Ast+Ap)/Ag

!   Calculate Cracking Moment (Mcr)

WorkVar1 = (1  + 50*r)
fcs = 1.5*r*Ep*us*1E-6/WorkVar1

Mcr = 1E-6*Zgb*(f`cf - fcs + 1E3*P/Ag) + P*e*1E-3  

!7b  SHEAR DESIGN - Flexure shear cracking

scpf = 1000*P/Ag + 1000*P*e*Yb/Ig

Mo = scpf*Zgb/10^6          ! Decompression moment
Vo = Mo*Vult/Mult

!   Calc dist from extreme comprn fibre to centroid of outermost layer of tensile R/F or tendons
!

SORT_LOW_HIGH (Ybarri, Fbari, Dsorti, Fsorti)
WorkVar1 = Dsort1

IF (do = 0) THEN do = D - WorkVar1

doMin = 0.8*D
Bw = 2*Tw

Beta1 = 1.1*(1.6 - do/1000)
IF (Beta1 < 1.1) THEN Beta1 = 1.1
Beta2 = 1.0     !  For Super-Ts (AS5100 Clause 8.2.7.1)

IF (Beta3 = 0) THEN Beta3 = 2*do*1000/x
IF (Beta3 < 1) THEN Beta3 = 1.0
IF (Beta3 > 2) THEN Beta3 = 2.0

WorkVar1 = (Ast+Ap)*f`cg/(Tw*do)
WorkVar2 = WorkVar1 ^ 0.333
V`uc = Beta1*Beta2*Beta3*Tw*do*WorkVar2/1000.
Vucs = V`uc + Vo

!7c  SHEAR DESIGN - Web shear cracking at centroid of composite section (direct solution)

sc = -P*1000/Ag + P*e*1000*(Yc - Yb)/Ic
Qlbwc = Qna/(Ic*Bw)

IF (sc < st) THEN
   tc = (st*st-st*sc)^0.5     ! Shear stress tau
ELSE
   tc = 0.
END IF

WorkVar1 = 0.5*sc
WorkVar2 = WorkVar1*WorkVar1
s1c = (WorkVar2 + tc*tc)^0.5 + WorkVar1
Vtc = 1E-3*tc/Qlbwc

Mtc = Mult*Vtc/Vult     ! Moment capacity corresponding to shear capacity Vtc

IF (Mtc > Mcr) THEN     !  Section is cracked if Mtc > cracking moment Mcr
    VuccNte$ = Cracked
    Vucc = 0
ELSE
    VuccNte$ = OK
    Vucc = Vtc
END IF

!7d  SHEAR DESIGN - Web shear cracking at the bottom of the web-flange interface
!
!  The algorithm to determine the ultimate concrete shear capacity,
!  Vtf, is based on RF Warner, BV Rangan, AS Hall and KA Faulkes,
!  "Concrete Structures" (Addison Wesley Longman, Aust).  Assume:
!
!   y2 = Dist from web/flange interface to composite centroid
!   y3 = Dist from web/flange interface to centroid of precast girder
!  Dwf = Distance from bottom of girder to web/flange interface
!    P = Final prestress force
!    e = Eccentricity (CG girder-strand group)
!   Ag = Area of precast girder
!   Ig = Moment of Inertia of precast girder
!   Ic = Composite moment of inertia
!   sf = Flexural stress at web-flange interface for uncracked section
!   st = Allowable principal tensile strength
!   tf = Ultimate shear stress capacity of concrete on its own (tau)
!   Bw = Sum of widths of both webs
!  Qfw = Shear flow constant at flange/web junction
!  Mcf = Corresponding moment factor (= Mult/Vult)
!        (used to calculate the corresponding flexural stress 
!         when the shear capacity is reached)
!
!  The solution for Vtf is based on the following relationships:
!
!  sf = (-P/Ag) - (P*e*y3/Ig) + Vtf*Mcf*y2/Ic .......... (1) Flexural stress
!
!  st = 0.33SQRT(f'cg) = SQRT((sf/2)^2 + tf^2) + sf/2 .. (2) Tensile stress capacity
!
!  Vtf = tf*Ic*Bw/Qfw                      ............. (3) Shear force capacity
!
!  Equations (1), (2) and (3) can be combined to form a quadratic with
!  the only unknown being Vtf.  The equations may be simplified by
!  grouping the constants  viz:
!
!  Equation (1):  Let
!    s1 = (-P*1000/Ag) - (P*e*y3*1000/Ig)
!    Mcf = Mult/Vult
!    s2 =  y2*10^6/Ic
!    s3 =  Mcf*s2
!        
!  Substituting, Equation (1) becomes:  sf = s1 + Vtf*s3
!
!  Equation (2):  Let  st = 0.33SQRT(f'cg)
!    then:               SQRT((sf/2)^2 + tf^2) + sf/2 = st
!    reorganise:         SQRT((sf/2)^2 + tf^2) = st - sf/2
!    square both sides:  (sf/2)^2 + tf^2 = st^2 + (sf/2)^2 - 2*st*sf/2
!    simplify:           tf^2 = st^2 - st*sf
!    reorganise:         tf^2 + st*sf - st^2 = 0.
!
!  Equation (3):  If  Qlbwf = Qfw/(Ic*Bw)  then  tf = Vtf*Qlbwf
!
!  Substitute for tf and sf in equation (2):
!               (Vtf*Qlbwf)^2 + st*(s1 + Vtf*s3) - st^2 = 0
!  Expand: Qlbwf^2 * Vtf^2 + st*s3*Vtf + st*s1 - st^2 = 0
!
!  To simplify the above expression let: 
!        a = Qlbwf^2
!        b = st*s3 = st*Mcf*s2
!        c = st*s1 - st*st
!
!  Then:   a*Vtf*Vtf + b*Vtf + c = 0  and the solution for Vtf is:
!             Vtf = (-b +- SQRT(b*b - 4*a*c))/2*a
!
!  Setting the determinant of this quadratic equation to:
!        f = b*b - 4*a*c
!  there is no solution if 'f' is negative.  If 'f' is positive then:
!      Vtf = (-b + SQRT(f))/2*a

y2 = Yc - Dwf
y3 = Yb - Dwf

Qlbwf = Qfw/(Ic*Bw)
Mcf = Mult/Vult

s1 = -P*1000/Ag - P*e*1000*y3/Ig    ! Convert kN to N
s2 = 10^6*y2/Ic                               ! Convert kN.m to N.mm

a = Qlbwf*Qlbwf*1E6
b = st*Mcf*s2
c = st*s1 - st*st

f = (b*b - 4*a*c)

IF (f < 0.) THEN
    s1f = 0.
    Vtf = 0.
    Mtf = 999999.
ELSE
    Vtf = (-b + f^0.5 )/(2*a)    !  Convert N to kN
    sf = -P*1000/Ag - P*e*1000*y3/Ig + Vtf*Mcf*y2*1E6/Ic
    tf = 1000*Vtf*Qlbwf
    WorkVar1 = (0.5*sf)
    WorkVar2 = WorkVar1*WorkVar1
    WorkVar3 = WorkVar2 + tf*tf
    WorkVar4 = WorkVar3^0.5
    s1f = WorkVar4 + 0.5*sf
    Mtf = Mult*Vtf/Vult       ! Moment capacity corresponding to shear capacity Vtf
END IF

!   Section is cracked if Mtf > cracking moment Mcr

IF (Mtf > Mcr) THEN
    Vucf = 0.
    VuccNte$ = Cracked
ELSE
    Vucf = Vtf
    VuccNte$ = OK
END IF


!   SHEAR DESIGN - Determine Vuc (minimum of Vucc, Vucs, Vucf)

Vuc = Vucs

IF (Vuc > Vucc) THEN Vuc = Vucc
IF (Vuc = 0)    THEN Vuc = Vucc
IF (Vuc > Vucf) THEN Vuc = Vucf
IF (Vuc = 0)    THEN Vuc = Vucf

!7e  SHEAR DESIGN - Design of web shear reinforcement

øVumin = 0.7*(Vuc + 0.6*Bw*do/1000)
øVumax = 0.7*0.2*f`cg*Bw*do/1000
OVumin = øVumin             ! HTML report can`t handle `ø` symbol
OVumax = øVumax            ! HTML report can`t handle `ø` symbol

Sv = Svmax
Svr1 = 0.
Svr2 = 0.
Asvsr = 0.
Asvsr1 = 0.
Asvsr2 = 0.

IF (øVumin > Vult) THEN
   Ov = 30.
   Svr = Svmax                     ! Required shear R/F spacing
   Vus = 0.
   AsNote1$ = Asv-min
   AsNote2$ = N/A
ELSE
   Ov = 0.
   Svr = 0.
   Vus = (Vult - 0.7*Vuc)/0.7
   AsNote1$ = N/A
   AsNote2$ = Asv Reqd
END IF

Ov1 = 30.
Ov2 = 30 + 15*(Vult - øVumin)/ (øVumax - øVumin)

IF (Ov2 > Ov1) THEN Ov1 = Ov2
IF (Ov1 > 45) THEN Ov1 = 45

IF (øVumin > Vult) THEN Ov1 = 0
IF (øVumin > Vult) THEN Ov2 = 0

IF (Dsr = 0) THEN Dsr = 16

Asv = 2*3.14*Dsr*Dsr/4

IF (fsyf = 0.) THEN fsyf = 250

IF (øVumin > Vult) THEN
    Ovc = Ov
    Asvsr1 = 0.35*Bw/fsyf     ! Required area of shear R/F per unit length
    Svr1 = Asv/Asvsr1
ELSE
    Ovc = Ov2
    WorkVar1 = Ov1*3.14/180         ! Convert to radians [for TAN function]
    WorkVar2 = TAN(WorkVar1)
    WorkVar3 =  ABS(WorkVar2)
    Asvsr2 = Vus*1000*WorkVar3/(do*fsyf)
    Svr2 = Asv/Asvsr2
END IF

Asvsr = Asvsr1
IF (Asvsr < Asvsr2) THEN Asvsr = Asvsr2
IF (Svr2 > 300) THEN Svr2 = Svmax

IF (øVumin > Vult) THEN
    Svr = Svmax
    Asvs = Asv/Sv
ELSE
    Svr = ABS(Svr2)
    Asvs = Asvsr2
END IF


!8a  DEFORMATIONS - Axial shortening and transfer hog

Xe = Pt*Lg*1E6/(Ag*Egmt)      ! Elastic chortening
u3 = k1s*850
Xs = u3*Lg/1000                   ! Shortening due to shrinkage
øccsc = 2*k2sc*k3sc
u4 = fcscgs*øccsc*1E6/Egmi
Xc = u4*Lg/1000                   ! Shortening due to creep

Xa = Xe + Xs + Xc                 ! Total axial shortening
Occsc = øccsc                        ! HTML reports can't handle ø symbol

Dsw = 40.*Msw * Lg * Lg * 1E12 / (384. * Egmt * Ig)
Hps = -Pt * 1000 * e * (Lg * 1000) * (Lg * 1000) / (8. * Egmt * Ig)
Htr = Dsw + Hps

!8b DEFORMATIONS - Deflections

Ddeck = 40. * Mslab * Lg*Lg * 1E12 / (384. * Egmi * Ig)
Dsdl  = 40. * Msdl  * Lg*Lg * 1E12 / (384. * Eg * Ic)
Ddl = Ddeck + Dsdl
Dratio = Lg*1000 / Dll
Dperm = Lg*1000 / Drperm

!8c DEFORMATIONS - Hog due to creep

WorkVar5 = Ti     !  Forces Ti to be printed in the output report
th1 = 2 * Ag / Per
Fratio3 = f`cmi / f`cg
øcc1 = øccb1 * k2c * k3c
Occ1 = øcc1                    ! HTML report can`t handle ø symbol
Occb1= øccb1                ! HTML report can`t handle ø symbol
Ecc1 = fcgs * øcc1 * 1E6 / Egmi
Pc1 = Ecc1 * Ep * Ap / 1E9

!8d DEFORMATIONS - Girder hog at installation

WorkVar5 = Ty     !  Forces Ty to be printed in the output report
Hgic = øcc1*Htr
Hgips = -(Pt + Pc1)*e*Lg*Lg*1E9/(8*Egmi*Ig)
Hgi = Hgic + Hgips + Dsw

!8e DEFORMATIONS - Final girder hog after installation

th2 = 2 * Ac / (Gp + 0.5 * Vp)
Fratio4 = f`cg / f`cg
øcc3 = øccb3 * (k2f-k2d) * k3f
Occ3 = øcc3                          ! HTML report can`t handle ø symbol
Occb3 = øccb3                     ! HTML report can`t handle ø symbol

Hel = Hgips + Dsw + Ddl
Hgfps = -P*e*Lg*Lg*1E9/(8*Egmt*Ig)
Hgfc = øcc3*Hel
Hgf = Hgfps + Hgfc + Dsw + Ddl