c///////////////////////////////////////////////////////////////////////
c
c  Virginia Tech Department of Aerospace and Ocean Engineeirng
c
c  Program LDstab estimates Lateral/Directional Stability and Control
c  derivatives for transonic tranport class aircraft. It consists of 
c  a main program stabin that reads in the data, and a subroutine stab
c  that does the calculations. 
c
c  The program was written by Joel Grasmeyer, and is described in
c  VPI-AOE-254, which explains the input data required. The sign
c  conventions are also described there.
c
c  program stabin
c
c    This program reads the stability input file, and calls the stab
c  subroutine to calculate the maximum available yawing moment
c  coefficient.
c
c  See stab.f for variable definitions
c
c  Created by:  Joel Grasmeyer
c  Last Modified:  10/28/97
c
c  Mods:
c	- Modified for use in Aircraft design class by W.H. Mason, 3/30/2000
c	- Changed header, added pause at end, added an output file,
c	  L.T.Leifsson March 2004
c	- Added Cnbeta contribution from wing using DATCOM method,
c       this required adding two input variables: xbar and cbar,
c	  by L.T. Leifsson, March 2004.
c
c///////////////////////////////////////////////////////////////////////

      program stabin

      implicit none

      character*72 infile, outfile, header, title
      integer i, write_flag, unit_in, unit_outfile, new, nef
      real dihedral_wing, z_wing, dia_fuse, sref, hspan_vtail,
     & depth_fuse_vtail, c_vtail_root, c_vtail_tip, mach_eo, l_vtail,
     & sweep_wing_1_2, hspan_wing, cl, z_vtail,
     & length_fuse, rho_eo, a_eo, mu_eo, cn_avail, dia_nacelle,
     & ar_vtail_eff, dr_max, sh, sweep_vtail_1_4, xbar, cbar

c Get input and output filenames from command line

c      write(6,10)
	write(6,3000)
	write(6,3001)
	write(6,3002)
	write(6,3003)
	write(6,3004)
	write(6,3005)
      read(5,20) infile
	write(6,3006)
	read(5,20) outfile

c Read input that is common to both modes
      unit_in = 10
      open(unit_in,file=infile)
      read(unit_in,"(a72)") header
      read(unit_in,"(a72)") title
      read(unit_in,*) write_flag
      read(unit_in,*) dihedral_wing
      read(unit_in,*) z_wing
      read(unit_in,*) dia_fuse
      read(unit_in,*) sref
      read(unit_in,*) hspan_vtail
      read(unit_in,*) depth_fuse_vtail
      read(unit_in,*) c_vtail_root
      read(unit_in,*) c_vtail_tip
      read(unit_in,*) mach_eo
      read(unit_in,*) sweep_vtail_1_4
      read(unit_in,*) sweep_wing_1_2
      read(unit_in,*) hspan_wing
      read(unit_in,*) cl
      read(unit_in,*) z_vtail
      read(unit_in,*) l_vtail
      read(unit_in,*) length_fuse
      read(unit_in,*) new
      read(unit_in,*) nef
      read(unit_in,*) dia_nacelle
      read(unit_in,*) sh
      read(unit_in,*) rho_eo
      read(unit_in,*) a_eo
      read(unit_in,*) mu_eo
      read(unit_in,*) dr_max
	read(unit_in,*) xbar
	read(unit_in,*) cbar

c Close input file
      close(unit_in)
c Open output file and write header
	unit_outfile = 20
	open(unit_outfile,file=outfile)
	write(unit_outfile,3000)
	write(unit_outfile,3001)
	write(unit_outfile,3002)
	write(unit_outfile,3003)
	write(unit_outfile,3004)

c Echo input to output file

c Call stab subroutine
      call stab(outfile,title,write_flag,dihedral_wing,z_wing,
     & dia_fuse,sref,hspan_vtail,depth_fuse_vtail,c_vtail_root,
     & c_vtail_tip,mach_eo,sweep_vtail_1_4,sweep_wing_1_2,hspan_wing,
     & cl,z_vtail,l_vtail,length_fuse,new,nef,dia_nacelle,
     & sh,rho_eo,a_eo,mu_eo,dr_max,xbar,cbar,ar_vtail_eff,cn_avail)

c Pause before closing program
	PAUSE 'Press ENTER to close program'
c Write statements and formats  
   10 format(/2x, 'Lateral-Directional Stability and Control Derivative',
     1            ' Estimation'/2x,
     1            '-intended for transonic tranports, i.e., B747 -'//2x,
     1            'Joel Grasmeyer'/2x,
     2            'Aerospace and Ocean Engineering, Virginia Tech'/2x,
     3            'Blacksburg, VA 24061'//2x,
     4            'Aircraft Design Class Software'//2x,
     5            'Enter the name of input data file')
   20 format(a20)
 3000 format('  **************************************************'/
     1        '    VT Aerospace Aircraft Design Software Series'/
     2		'  **************************************************'/)
 3001 format('                   Program LDstab'/
     1		'  The program estimates Lateral/Directional Stability'/ 
     2		'  and Control derivatives for transonic tranport class'/
     3		'  aircraft.'/)
 3002 format('  Written by Joel Grasmeyer.'/
     1		'  This program modified by'/
	2		'    - W.H. Mason, March 2000,'/
     3		'    - L.T. Leifsson, March 2004.'/)
 3003 format('  The Department of Aerospace and Ocean Engineering,'/
     1		'  Virginia Tech, Blacksburg, VA 24061'/
     2		'  http://www.aoe.vt.edu, email: whmason@vt.edu'/)
 3004 format('  **************************************************'/)
 3005	format('  Enter input filename:')
 3006	format('  Enter output filename:')
      end ! End program stabin

c///////////////////////////////////////////////////////////////////////
c
c  subroutine stab
c
c    This subroutine calculates the maximum available yawing moment
c  coefficient.  Note that right rudder deflection is defined as
c  positive, and right aileron up, left aileron down is defined as
c  positive.  Both of these control deflections generate positive
c  moments about their respective axes.  This is the convention used
c  by Roskam.
c
c  Inputs
c
c  outfile      output filename
c  title        title of aircraft configuration
c  write_flag   write flag (0 = no output file, 1 = output file written)
c  dihedral_wing  wing dihedral angle (deg)
c  z_wing       distance from body centerline to quarter-chord point of
c                exposed wing root chord, positive for the quarter-chord
c                point below the body centerline (ft)
c  dia_fuse     fuselage diameter (ft)
c  sref         wing reference area (ft^2)
c  hspan_vtail  vertical tail span (ft)
c  depth_fuse_vtail  fuselage depth at the fuselage station of the
c                     quarter-chord of the vertical tail (ft)
c  c_vtail_root root chord of vertical tail
c  c_vtail_tip  tip chord of vertical tail
c  mach_eo      mach number
c  sweep_vtail_1_4  vertical tail quarter-chord sweep angle (deg)
c  sweep_wing_1_2  average wing half-chord sweep angle (deg)
c  hspan_wing   wing half-span (ft)
c  cl           lift coefficient
c  z_vtail      vertical distance from CG to AC of vertical tail (ft)
c  l_vtail      horizontal distance from CG to AC of vertical tail (ft)
c  length_fuse  fuselage length (ft)
c  new          number of engines on the wing
c  nef          number of engines on the fuselage
c  dia_nacelle  nacelle diameter (ft)
c  rho_eo       density at engine-out flight condition (slug/ft^3)
c  a_eo         speed of sound at engine-out flight condition (ft/s)
c  mu_eo        viscosity at engine-out flight condition (slug/(ft-s))
c  dr_max       maximum allowable steady-state rudder deflection (deg)
c  xbar		  longitudinal distance (positive rearward) from coordinate
c			  origin (usually the center of gravity) to the wing aero-
c			  dynamic center (ft), L.T. Leifsson March 2004.
c  cbar		  mean aerodynamic chord (ft), L.T. Leifsson March 2004.
c
c  Outputs
c
c  ar_vtail_eff  effective aspect ratio of vertical tail
c  cn_avail     maximum available yawing moment coefficient
c
c  Internal Variables
c
c  alpha        angle of attack (rad)
c  alpha_d      section lift effectiveness
c  alpha_d_cl   section flap effectiveness (from Figure 10.2)
c  ar           wing aspect ratio
c  ar_vtail     actual aspect ratio of vertical tail
c  beta         sideslip angle, positive from the right (rad)
c  beta_m       square root of (1 - mach_eo)**2
c  cf_c         ratio of flap chord to wing or tail chord
c  cl_alpha     lift-curve slope of entire aircraft (rad^-1)
c  cl_alpha_2d  2-dimensional lift-curve slope at MAC (rad^-1)
c  cl_alpha_vtail  original lift-curve slope of vertical tail (rad^-1)
c  cl_alpha_vtail_eff  effective lift-curve slope of vertical tail (rad^-1)
c  cl_beta      variation of rolling moment coefficient with sideslip angle
c  cl_beta_htail  horizontal tail contribution to cl_beta
c  cl_beta_vtail  vertical tail contribution to cl_beta
c  cl_beta_wingbody  wing-body contribution to cl_beta
c  cl_d         rolling effectiveness of partial-chord controls
c  cl_da        variation of rolling moment coefficient with aileron deflection
c  cl_dr        variation of rolling moment coefficient with rudder deflection
c  cn_beta      variation of yawing moment coefficient with sideslip angle
c  cn_beta_fuse  fuselage contribution to cn_beta
c  cn_beta_vtail  vertical tail contribution to cn_beta
c  cn_beta_wing  wing contribution to cn_beta
c  cn_da        variation of yawing moment coefficient with aileron deflection
c  cn_dr        variation of yawing moment coefficient with rudder deflection
c  cy_beta      variation of side force coefficient with sideslip angle
c  cy_beta_fuse  fuselage contribution to cy_beta
c  cy_beta_vtail  vertical tail contribution to cy_beta
c  cy_beta_wing  wing contribution to cy_beta
c  cy_da        variation of side force coefficient with aileron deflection
c  cy_dr        variation of side force coefficient with rudder deflection
c  da           aileron deflection, positive for right aileron up, left
c                aileron down (rad)
c  dr           rudder deflection, positive for right deflection (rad)
c  eqn_7_5      factor estimated by Equation 7.5
c  eqn_7_10     body-induced effect on wing height (from Equation 7.10)
c  eqn_7_11     d in Equation 7.10 (estimated from Equation 7.11)
c  eqn_7_12     another body-induced effect on wing height (from Equation 7.12)
c  eqn_11_2     rolling effectiveness of two full-chord ailerons (Eqn. 11.2)
c  f_cy_beta    correction factor for cy_beta
c  f_cl_beta    correction factor for cl_beta
c  f_cn_beta    correction factor for cn_beta
c  f_cl_da      correction factor for cl_da
c  f_cn_da      correction factor for cn_da
c  f_cy_dr      correction factor for cy_dr
c  f_cl_dr      correction factor for cl_dr
c  f_cn_dr      correction factor for cn_dr
c  fig_7_1      wing-body interference factor from Figure 7.1
c  fig_7_3      empirical factor from Figure 7.3
c  fig_7_5      ratio of the aspect ratio of the vertical panel in the presence
c                of the body to that of the isolated panel (from Figure 7.5)
c  fig_7_6      ratio of the vertical panel aspect ratio in the presence of
c                the horizontal tail and body to that of the panel in the
c                presence of the body alone (from Figure 7.6)
c  fig_7_7      factor accounting for relative size of horizontal and vertical
c                tails (from Figure 7.7)
c  fig_7_11     wing sweep contribution to cl_beta_wingbody (from Figure 7.11)
c  fig_7_12     compressibility correction to wing sweep (from Figure 7.12)
c  fig_7_13     fuselage correction factor (from Figure 7.13)
c  fig_7_14     aspect ratio contribution to cl_beta_wingbody (from Figure 7.14)
c  fig_7_15     dihedral effect on cl_beta (from Figure 7.15)
c  fig_7_16     compressibility correction to dihedral effect (from Figure 7.16)
c  fig_7_19     factor for body and body + wing effects (from Figure 7.19)
c  fig_7_20     Reynold's number factor for the fuselage (from Figure 7.20)
c  fig_10_2     flap chord factor (from Figure 10.2)
c  fig_10_3     span factor for plain flap (from Figure 10.3)
c  fig_10_5     theoretical lift effectiveness of plain TE flaps (Fig. 10.5)
c  fig_10_6     empirical correction for plain TE flaps (Fig. 10.6)
c  fig_10_7     empirical correction for lift effectiveness of plain flaps
c                at high flap deflections (from Figure 10.7)
c  fig_11_1     rolling moment effectiveness parameter (from Figure 11.1)
c  flap_eff_ratio  flap effectiveness ratio (from Figure 10.2)
c  i            index
c  k            empirical factor for estimating the variation of yawing
c                moment coefficient with aileron deflection
c  kappa        ratio of the actual lift-curve slope to 2*pi
c  phi          bank angle, positive to the right (rad)
c  q            dynamic pressure (lb/ft^2)
c  re_fuse      fuselage Reynolds number
c  sbs          body side area (ft^2)
c  sh           area of horizontal tail (ft^2)
c  sv           area of vertical tail (ft^2)
c  sweep_vtail_1_2  vertical tail half-chord sweep angle (deg)
c  x            temporary variable for curve fits
c
c  Created by:  Joel Grasmeyer
c  Last Modified:  10/28/97
c
c  I/O changes by W.H. Mason, 3/30/2000
c
c///////////////////////////////////////////////////////////////////////

      subroutine stab(outfile,title,write_flag,dihedral_wing,z_wing,
     & dia_fuse,sref,hspan_vtail,depth_fuse_vtail,c_vtail_root,
     & c_vtail_tip,mach_eo,sweep_vtail_1_4,sweep_wing_1_2,hspan_wing,
     & cl,z_vtail,l_vtail,length_fuse,new,nef,dia_nacelle,
     & sh,rho_eo,a_eo,mu_eo,dr_max,xbar,cbar,ar_vtail_eff,cn_avail)

      implicit none

      character*72 outfile, title
      integer i, write_flag, unit_out, unit_outfile, new, nef
      real pi, dihedral_wing, z_wing, dia_fuse, sref, hspan_vtail, ar,
     & depth_fuse_vtail, c_vtail_root, c_vtail_tip, mach_eo, sv, sh,
     & sweep_wing_1_2, hspan_wing, cy_beta, ar_vtail, k,
     & cy_beta_wing, cy_beta_fuse, cy_beta_vtail, ar_vtail_eff, alpha,
     & cl_alpha_vtail, beta_m, eqn_7_5, kappa, cl_alpha_vtail_eff, q,
     & cl_beta, cl_beta_htail, cl_beta_vtail, cl_beta_wingbody, sbs,
     & cl, fig_7_1, fig_7_5, fig_7_6, fig_7_7, fig_7_11, cn_da,
     & fig_7_12, fig_7_13, fig_7_14, fig_7_15, fig_7_16, eqn_7_10, da,
     & eqn_7_11, eqn_7_12, cl_alpha, z_vtail, l_vtail, cn_beta_fuse,
     & cn_beta_vtail, cn_beta, cn_beta_wing, fig_7_19, fig_7_20, phi,
     & length_fuse, re_fuse, cl_da, cy_da, fig_11_1, eqn_11_2, cl_d,
     & alpha_d, fig_10_5, fig_10_6, cl_alpha_2d, cl_dr, dr,
     & cn_dr, cy_dr, fig_10_2, fig_10_7, alpha_d_cl, flap_eff_ratio,
     & fig_10_3, cf_c, beta, rho_eo, a_eo, mu_eo, cn_avail, fig_7_3,
     & dia_nacelle, dr_max, f_cy_beta, f_cl_beta,
     & f_cn_beta, f_cl_da, f_cn_da, f_cy_dr, f_cl_dr, f_cn_dr, x,
     & sweep_vtail_1_4, sweep_vtail_1_2,
     & cnbw1, cnbw2, cnbw3, cnbw4, cnbw5, cnbw6, cn_beta_wingterm, 
     & cbar, xbar

      pi = acos(-1.)

c Convert sweep angles from degrees to radians
cwhm THIS IS A VERY DANGEROUS PRACTICE!! DON'T USE AS AN EXAMPLE
cwhm of acceptable coding style
      sweep_wing_1_2 = sweep_wing_1_2*pi/180.
      sweep_vtail_1_4 = sweep_vtail_1_4*pi/180.

c Write input data to screen for confirmation
      if (write_flag .eq. 1) then
        unit_out = 6
        write(unit_out,90) title
        write(unit_out,*)
        write(unit_out,"(a25)") '   Input Data'
        write(unit_out,*)
        write(unit_out,101) write_flag, '= write_flag'
        write(unit_out,100) dihedral_wing, '= dihedral_wing'
        write(unit_out,100) z_wing, '= z_wing'
        write(unit_out,100) dia_fuse, '= dia_fuse'
        write(unit_out,100) sref, '= sref'
        write(unit_out,100) hspan_vtail, '= hspan_vtail'
        write(unit_out,100) depth_fuse_vtail, '= depth_fuse_vtail'
        write(unit_out,100) c_vtail_root, '= c_vtail_root'
        write(unit_out,100) c_vtail_tip, '= c_vtail_tip'
        write(unit_out,100) mach_eo, '= mach_eo'
        write(unit_out,100) sweep_vtail_1_4*180./pi, '= sweep_vtail_1_4'
        write(unit_out,100) sweep_wing_1_2*180./pi, '= sweep_wing_1_2'
        write(unit_out,100) hspan_wing, '= hspan_wing'
        write(unit_out,100) cl, '= cl'
        write(unit_out,100) z_vtail, '= z_vtail'
        write(unit_out,100) l_vtail, '= l_vtail'
        write(unit_out,100) length_fuse, '= length_fuse'
        write(unit_out,101) new, '= new'
        write(unit_out,101) nef, '= nef'
        write(unit_out,100) dia_nacelle, '= dia_nacelle'
        write(unit_out,100) sh, '= sh'
        write(unit_out,100) rho_eo, '= rho_eo'
        write(unit_out,100) a_eo, '= a_eo'
        write(unit_out,103) mu_eo, '= mu_eo'
        write(unit_out,100) dr_max, '= dr_max'
	  write(unit_out,100) xbar, '= xbar'
	  write(unit_out,100) cbar, '= cbar'
      end if

c Echo input to output file
      unit_outfile = 20
	write(unit_outfile,90) title
      write(unit_outfile,*)
      write(unit_outfile,"(a25)") '   Input Data'
      write(unit_outfile,*)
      write(unit_outfile,101) write_flag, '= write_flag'
      write(unit_outfile,100) dihedral_wing, '= dihedral_wing'
      write(unit_outfile,100) z_wing, '= z_wing'
      write(unit_outfile,100) dia_fuse, '= dia_fuse'
      write(unit_outfile,100) sref, '= sref'
      write(unit_outfile,100) hspan_vtail, '= hspan_vtail'
      write(unit_outfile,100) depth_fuse_vtail, '= depth_fuse_vtail'
      write(unit_outfile,100) c_vtail_root, '= c_vtail_root'
      write(unit_outfile,100) c_vtail_tip, '= c_vtail_tip'
      write(unit_outfile,100) mach_eo, '= mach_eo'
      write(unit_outfile,100) sweep_vtail_1_4*180./pi,'=sweep_vtail_1_4'
      write(unit_outfile,100) sweep_wing_1_2*180./pi, '= sweep_wing_1_2'
      write(unit_outfile,100) hspan_wing, '= hspan_wing'
      write(unit_outfile,100) cl, '= cl'
      write(unit_outfile,100) z_vtail, '= z_vtail'
      write(unit_outfile,100) l_vtail, '= l_vtail'
      write(unit_outfile,100) length_fuse, '= length_fuse'
      write(unit_outfile,101) new, '= new'
      write(unit_outfile,101) nef, '= nef'
      write(unit_outfile,100) dia_nacelle, '= dia_nacelle'
      write(unit_outfile,100) sh, '= sh'
      write(unit_outfile,100) rho_eo, '= rho_eo'
      write(unit_outfile,100) a_eo, '= a_eo'
      write(unit_outfile,103) mu_eo, '= mu_eo'
      write(unit_outfile,100) dr_max, '= dr_max'
	write(unit_outfile,100) xbar, '= xbar'
	write(unit_outfile,100) cbar, '= cbar'

c Calculate stability and control derivatives via Roskam's methods

c Sideslip angle derivatives
      cy_beta_wing = -0.0001*abs(dihedral_wing)*180./pi

c Estimate fig_7_1 from Figure 7.1 (curve fit)
      if (z_wing/(dia_fuse/2.) .le. 0.) then
        fig_7_1 = 0.85*(-z_wing/(dia_fuse/2.)) + 1.
      elseif (z_wing/(dia_fuse/2.) .gt. 0.) then
        fig_7_1 = 0.5*z_wing/(dia_fuse/2.) + 1.
      end if

c Estimate the side force coefficient due to the fuselage and nacelles
      cy_beta_fuse = -2.*fig_7_1*( pi*(dia_fuse/2.)**2 +
     & (new + nef)*pi*(dia_nacelle/2.)**2 )/sref

c Estimate fig_7_3 from Figure 7.3 (curve fit)
      x = hspan_vtail/depth_fuse_vtail
      if (x .le. 2.) then
        fig_7_3 = 0.75
      elseif (x .gt. 2. .and. x .lt. 3.5) then
        fig_7_3 = x/6. + 5./12.
      elseif (x .ge. 3.5) then
        fig_7_3 = 1.
      end if

c Estimate fig_7_5 from Figure 7.5 (curve fit for taper ratio <= 0.6)
      x = hspan_vtail/depth_fuse_vtail
      fig_7_5 = 0.002*x**5 - 0.0464*x**4 + 0.404*x**3 - 1.6217*x**2 +
     & 2.7519*x + 0.0408

c Factor from Figure 7.6 is for zh/bv = 0.
      fig_7_6 = 1.1

c Estimate fig_7_7 from Figure 7.7 (curve fit)
      sv = hspan_vtail*(c_vtail_root + c_vtail_tip)/2.
      x = sh/sv
      fig_7_7 = -0.0328*x**4 + 0.2885*x**3 - 0.9888*x**2 + 1.6554*x -
     & 0.0067

c Estimate the effective aspect ratio for the vertical tail
      ar_vtail = hspan_vtail**2/sv
      ar_vtail_eff = fig_7_5*ar_vtail*(1. + fig_7_7*(fig_7_6 - 1.))

c Assume the section lift-curve slope is 2.*pi
      cl_alpha_vtail = 2.*pi

c Estimate the effective lift-curve slope for the vertical tail
      kappa = cl_alpha_vtail/(2.*pi)
      beta_m = sqrt( 1. - mach_eo**2 )
      sweep_vtail_1_2 = atan( (c_vtail_root/4. + hspan_vtail*
     & tan(sweep_vtail_1_4) + c_vtail_tip/4. - c_vtail_root/2.)/
     & hspan_vtail )
      cl_alpha_vtail_eff = 2*pi*ar_vtail_eff/( 2. +
     & sqrt( ar_vtail_eff**2*beta_m**2/kappa**2*( 1. +
     & tan(sweep_vtail_1_2)**2/beta_m**2 ) + 4. ) )

c Estimate the third term in eqn. 7.4 from eqn. 7.5
      eqn_7_5 = 0.724 + 3.06*sv/sref/(1. + cos(sweep_vtail_1_4)) +
     & 0.4*z_wing/dia_fuse + 0.009*ar_vtail_eff
      cy_beta_vtail = -fig_7_3*cl_alpha_vtail_eff*eqn_7_5*sv/sref

c Calculate total variation of side force coefficient with sideslip angle
      cy_beta = cy_beta_wing + cy_beta_fuse + cy_beta_vtail

c Factor from Figure 7.11 is approximated by a curve fit for lambda = 0.5
      fig_7_11 = -0.004/45*sweep_wing_1_2*180./pi

c Factor from Figure 7.12 is approximated for 747 and 777 configurations
c at low Mach numbers
      fig_7_12 = 1.0

c Factor from Figure 7.13 is approximated for 747 and 777 configurations
      fig_7_13 = 0.85

c Factor from Figure 7.14 is approximated for lambda = 0.5 and high AR
      fig_7_14 = 0.000

c Factor from Figure 7.15 is approximated by a linear curve fit for
c  lambda = 0.5, low sweep, and high AR
      ar = (2.*hspan_wing)**2/sref
      fig_7_15 = -0.00012 - 0.00013/10*ar

c Factor from Figure 7.16 is approximated for 747 and 777 configurations
c at low Mach numbers
      fig_7_16 = 1.0

c Estimate body-induced effect on wing height from eqns. 7.10, 7.11, and 7.12
      eqn_7_11 = sqrt(pi*(dia_fuse/2.)**2/0.7854)
      eqn_7_10 = -0.0005*sqrt(ar)*(eqn_7_11/(2.*hspan_wing))**2
      eqn_7_12 = -1.2*sqrt(ar)/(180./pi)*z_wing/(2.*hspan_wing)*
     & 2.*eqn_7_11/(2.*hspan_wing)

c Wing-body contribution to cl_beta (wing twist effect is neglected)
      cl_beta_wingbody = ( cl*(fig_7_11*fig_7_12*fig_7_13 +
     & fig_7_14) + dihedral_wing*(fig_7_15*fig_7_16 + eqn_7_10) +
     & eqn_7_12 )*180./pi

c Since the horizontal tail has a small lift coefficient, small dihedral,
c  and small area relative to the wing, it is negligible.
      cl_beta_htail = 0.

c Estimate the angle of attack from the CL by assuming cl_alpha = 5.7
c based on the value given for the 747 in Nelson.  The angle of attack
c of the fuselage reference line is given as 5.7 deg for the M=0.25
c flight condition in NASA CR-2144.  Therefore, the effective wing
c incidence angle (with flaps at 20 deg) is 5.457 deg.
      cl_alpha = 5.7
      alpha = cl/cl_alpha - 5.457*pi/180.

c Estimate the vertical tail contribution to cl_beta
      cl_beta_vtail = cy_beta_vtail*( z_vtail*cos(alpha) - l_vtail*
     & sin(alpha) )/(2.*hspan_wing)

c Calculate total variation of rolling moment coefficient with sideslip angle
      cl_beta = cl_beta_wingbody + cl_beta_htail + cl_beta_vtail

c Wing contribution to cn_beta is negligible for small angles of attack.
c      cn_beta_wing = 0.
c Calculate contribution of wing to cn_beta using DATCOM method.
c   Code added by L.T. Leifsson, code written by W.H. Mason, March 2004. 
	cnbw1 = 1.0/(4.0*pi*ar)
	cnbw2 = tan(sweep_wing_1_2)/pi/ar/(ar + 4.0*cos(sweep_wing_1_2))
	cnbw3 = cos(sweep_wing_1_2)
	cnbw4 = ar/2.0
	cnbw5 = ar**2.0/8.0/cos(sweep_wing_1_2)
	cnbw6 = 6.0*xbar/cbar*sin(sweep_wing_1_2)/ar
	cn_beta_wingterm = cnbw1 - cnbw2*(cnbw3 - cnbw4 - cnbw5 + cnbw6)
	cn_beta_wing = cl**2.0*cn_beta_wingterm

c Estimate empirical factor for body and body + wing effects from Figure 7.19
c  Constant value assumed for 747 and 777-like configurations
      fig_7_19 = 0.0011

c Calculate fuselage Reynolds number at the engine-out flight condition
      re_fuse = rho_eo*mach_eo*a_eo*length_fuse/mu_eo

c Estimate fuselage Reynolds number effect on wing-body from Figure 7.20
      fig_7_20 = 1. + 1.2/log(350.)*log(re_fuse/1000000.)

c Estimate fuselage contribution to cn_beta
      sbs = 0.83*dia_fuse*length_fuse
      cn_beta_fuse = -180./pi*fig_7_19*fig_7_20*sbs/sref*
     & length_fuse/(2.*hspan_wing)

c Estimate vertical tail contribution to cn_beta
      cn_beta_vtail = -cy_beta_vtail*( l_vtail*cos(alpha) +
     & z_vtail*sin(alpha) )/(2.*hspan_wing)

c Calculate total variation of rolling moment coefficient with sideslip angle
      cn_beta = cn_beta_wing + cn_beta_fuse + cn_beta_vtail

c Assume variation of sideforce coefficient with aileron deflection is zero
      cy_da = 0.

c Estimate the rolling moment effectiveness parameter from Figure 11.1
c  for lambda = 0.5, and for 747 and 777-like ailerons at mach 0.2
      fig_11_1 = 0.18

c Assume the 2D lift-curve slope is 2*pi/beta_m
      cl_alpha_2d = 2*pi/beta_m
      kappa = cl_alpha_2d/(2.*pi/beta_m)

c Estimate the rolling effectiveness of two full-chord controls by Eqn. 11.2
      eqn_11_2 = kappa/beta_m*fig_11_1

c Estimate aileron effectiveness by assuming cf/c = 0.20 and t/c = 0.08
      fig_10_5 = 3.5
      fig_10_6 = 1.0
      cl_d = fig_10_6*fig_10_5
      alpha_d = cl_d/cl_alpha_2d

c Determine the rolling effectiveness of the partial-chord controls by
c Eqn. 11.3.  Note that this is the change in cl with respect to a change
c in the sum of the left and right aileron deflections (d).
      cl_d = alpha_d*eqn_11_2

c Estimate variation of rolling moment coefficient with aileron deflection
c by neglecting differential control effects.  Since the aileron deflection
c (da) is defined as half of the sum of the left and right deflections, cl_d
c from the equation above must be divided by 2.
      cl_da = cl_d/2.

c The method in Roskam for estimating cn_da does not account for the
c effect of differential ailerons and the use of spoilers for roll control
c on the yaw moment.  Therefore, the factor k is estimated
c based on the ratio of cn_da to cl_da from the 747 flight test data
c presented in Nelson.  Note that the effect of cl is absorbed into
c the factor k.
      k = 0.0064/0.0461

c Estimate variation of yawing moment coefficient with aileron deflection
      cn_da = k*cl_da

c Estimate the flap chord factor from Figure 10.2 for cf/c = 0.33
c  The flap effectiveness ratio is estimated with a piecewise curve fit
      cf_c = 0.33
      alpha_d_cl = -sqrt( 1. - (1. - cf_c)**2 )
      if (alpha_d_cl .ge. -0.5) then
        flap_eff_ratio = 1.42 + 1.8*alpha_d_cl
      elseif (alpha_d_cl .ge. -0.6) then
        flap_eff_ratio = 1.32 + 1.6*alpha_d_cl
      elseif (alpha_d_cl .ge. -0.7) then
        flap_eff_ratio = 1.08 + 1.2*alpha_d_cl
      else
        flap_eff_ratio = 0.94 + alpha_d_cl
      end if
      flap_eff_ratio = 1. + flap_eff_ratio/( ar_vtail_eff -
     & 0.5*(-alpha_d_cl - 2.1) )
      fig_10_2 = flap_eff_ratio*alpha_d_cl

c Estimate empirical correction for lift effectiveness of plan flaps at
c from Figure 10.7 for cf/c = 0.33.
      x = dr_max
      if (x .lt. 15.) then
        fig_10_7 = 1.
      else
        fig_10_7 = 4e-7*x**4 - 7e-5*x**3 + 0.0047*x**2 - 0.1453*x +
     & 2.3167
      end if

c Estimate span factor for plain flap from Figure 10.3 for delta eta = 0.85
      fig_10_3 = 0.95

c Estimate variation of sideforce coefficient with rudder deflection
      cy_dr = cl_alpha_vtail_eff*fig_10_2*fig_10_7*fig_10_3*sv/sref

c Estimate variation of rolling moment coefficient with rudder deflection
      cl_dr = cy_dr*( z_vtail*cos(alpha) - l_vtail*sin(alpha) )/
     & (2.*hspan_wing)

c Estimate variation of yawing moment coefficient with rudder deflection
      cn_dr = -cy_dr*( l_vtail*cos(alpha) + z_vtail*sin(alpha) )/
     & (2.*hspan_wing)

c Multiply empirical estimates by their respective correction factors
c The correction factors are the ratio of the actual 747 derivatives to
c the 747 derivatives predicted by the method above at the M=0.25 flight
c condition defined in NASA CR-2144 and Nelson.  The rudder deflection
c was 15 deg for this calibration.
      cy_beta = 1.4068*cy_beta
      cl_beta = 0.7253*cl_beta
      cn_beta = 1.4232*cn_beta	! New correction factor, Cnbetawing now 
      cl_da   = 0.9202*cl_da		! estimated. Assuming xbar = 0.0 ft and
      cn_da   = 0.9143*cn_da		! cbar = 27.31 ft, for 747,Leifsson, March '04.
      cy_dr   = 0.6132*cy_dr
      cl_dr   = 0.3004*cl_dr
      cn_dr   = 0.7320*cn_dr

c Set the rudder deflection to 20 deg, and the bank angle to 5 deg
      dr = dr_max*pi/180.
      phi = 5.*pi/180.

c Solve for the sideslip angle and aileron deflection
      beta = -( cy_dr*dr + cl*sin(phi) )/cy_beta
      da = -( cl_dr*dr + cl_beta*beta)/cl_da

c Calculate the maximum available yawing moment coefficient
      cn_avail = cn_da*da + cn_dr*dr + cn_beta*beta

c Write output data
      if (write_flag .eq. 1) then
        write(unit_out,*)
        write(unit_out,"(a26)") '   Output Results'
        write(unit_out,96)
        write(unit_out,104) cy_beta_wing, '= cy_beta_wing'
        write(unit_out,104) cy_beta_fuse, '= cy_beta_fuse'
        write(unit_out,104) cy_beta_vtail, '= cy_beta_vtail'
        write(unit_out,*)
        write(unit_out,105) cy_beta, '= cy_beta'
        write(unit_out,*)
        write(unit_out,104) cl_beta_wingbody, '= cl_beta_wingbody'
        write(unit_out,104) cl_beta_htail, '= cl_beta_htail'
        write(unit_out,104) cl_beta_vtail, '= cl_beta_vtail'
        write(unit_out,*)
        write(unit_out,105) cl_beta, '= cl_beta'
        write(unit_out,*)
        write(unit_out,104) cn_beta_wing, '= cn_beta_wing'
        write(unit_out,104) cn_beta_fuse, '= cn_beta_fuse'
        write(unit_out,104) cn_beta_vtail, '= cn_beta_vtail'
        write(unit_out,*)
        write(unit_out,105) cn_beta, '= cn_beta'
        write(unit_out,97)
        write(unit_out,100) cy_da, '= cy_da'
        write(unit_out,100) cl_da, '= cl_da'
        write(unit_out,100) cn_da, '= cn_da'
        write(unit_out,*)
        write(unit_out,100) cy_dr, '= cy_dr'
        write(unit_out,100) cl_dr, '= cl_dr'
        write(unit_out,100) cn_dr, '= cn_dr'
        write(unit_out,98)
        write(unit_out,100) beta*180./pi, '= beta (deg)'
        write(unit_out,100) phi*180./pi, '= phi (deg)'
        write(unit_out,100) da*180./pi, '= da (deg)'
        write(unit_out,100) dr*180./pi, '= dr (deg)'
        write(unit_out,100) ar_vtail_eff, '= ar_vtail_eff'
        write(unit_out,100) cn_avail, '= cn_avail'
        write(unit_out,99)
      endif

c Write results to output file
      write(unit_outfile,*)
      write(unit_outfile,"(a26)") '   Output Results'
      write(unit_outfile,96)
      write(unit_outfile,104) cy_beta_wing, '= cy_beta_wing'
      write(unit_outfile,104) cy_beta_fuse, '= cy_beta_fuse'
      write(unit_outfile,104) cy_beta_vtail, '= cy_beta_vtail'
      write(unit_outfile,*)
      write(unit_outfile,105) cy_beta, '= cy_beta'
      write(unit_outfile,*)
      write(unit_outfile,104) cl_beta_wingbody, '= cl_beta_wingbody'
      write(unit_outfile,104) cl_beta_htail, '= cl_beta_htail'
      write(unit_outfile,104) cl_beta_vtail, '= cl_beta_vtail'
      write(unit_outfile,*)
      write(unit_outfile,105) cl_beta, '= cl_beta'
      write(unit_outfile,*)
      write(unit_outfile,104) cn_beta_wing, '= cn_beta_wing'
      write(unit_outfile,104) cn_beta_fuse, '= cn_beta_fuse'
      write(unit_outfile,104) cn_beta_vtail, '= cn_beta_vtail'
      write(unit_outfile,*)
      write(unit_outfile,105) cn_beta, '= cn_beta'
      write(unit_outfile,97)
      write(unit_outfile,100) cy_da, '= cy_da'
      write(unit_outfile,100) cl_da, '= cl_da'
      write(unit_outfile,100) cn_da, '= cn_da'
      write(unit_outfile,*)
      write(unit_outfile,100) cy_dr, '= cy_dr'
      write(unit_outfile,100) cl_dr, '= cl_dr'
      write(unit_outfile,100) cn_dr, '= cn_dr'
      write(unit_outfile,98)
      write(unit_outfile,100) beta*180./pi, '= beta (deg)'
      write(unit_outfile,100) phi*180./pi, '= phi (deg)'
      write(unit_outfile,100) da*180./pi, '= da (deg)'
      write(unit_outfile,100) dr*180./pi, '= dr (deg)'
      write(unit_outfile,100) ar_vtail_eff, '= ar_vtail_eff'
      write(unit_outfile,100) cn_avail, '= cn_avail'
      write(unit_outfile,99)

c Close output file
	close(unit_outfile)
c Write statements and formats
   90 format(/4x,'Data set title: ',a72)
   96 format(/3x,'stability derivatives'/)
   97 format(/3x,'control derivatives'/)
   98 format(/3x,'engine out information - 5 deg bank and max rudder'/)
   99 format(/3x,'Calculations completed')
  100 format(3x,f11.4, 1x, a)
  101 format(10x, i4, 1x, a)
  102 format(3x,f11.0, 1x, a)
  103 format(3x,g11.4, 1x, a)
  104 format(22x,f11.4, 1x, a)
  105 format(3x,f11.4, 1x, a, ' (adjusted)')

c Convert sweep angles from radians back to degrees
cwhm THIS IS A VERY DANGEROUS PRACTICE!! DON'T USE AS AN EXAMPLE
cwhm of acceptable coding style
      sweep_wing_1_2 = sweep_wing_1_2*180./pi
      sweep_vtail_1_4 = sweep_vtail_1_4*180./pi

      return
      end

