program loads2d
c
c*********************************************************************c
c*                                                                   *c
c*  Written by Tim J. Cowan                                          *c
c*                                                                   *c
c*********************************************************************c
c*
c*  Description: Uses trapezoidal integration to compute the 2D 
c*               aerodynamic loads for a cross section of Cp data.
c*
c*  References:  Bertin & Smith.  Aerodynamics for Engineers.  2nd 
c*               Edition.  Prentice Hall, Inc..  New Jersey: 1989.
c*
c*  Updated:     1/6/02 - First version.
c*
c*********************************************************************c
c*
c*  Input :   User provided data file containing x-coordinates, point
c*            slopes and measured pressure coefficients at stations
c*            along the surface of an airfoil.  If the points slope
c*            is unknown or undefined, set it to zero.  The point
c*            slope is required for an accurate form drag calculation.
c*            Header of input data file defines the free-stream
c*            angle of attack (alp) and the number of stations along
c*            the upper and lower surfaces (npu & npl). 
c*
c*  Output :  To screen, x & y force coefficients plus lift and form
c*            drag coefficients based on user defined free-stream
c*            angle of attack.
c*
c*********************************************************************c
c*
c*  Variable definitions
c*
c*  name  => description (default value) {additional information}
c*
c*  alp   => free-stream angle of attack in degrees
c*  npu   => number of points (stations) along the upper surface
c*  npl   => number of points (stations) along the lower surface
c*  nnd   => total number of points (stations) along the surface
c*
c*  X     => vector of x-coordinates for each station
c*  YP    => vector of point slopes for each station
c*  CP    => vector of pressure coefficients for each station
c*
c*********************************************************************c
c
      parameter ( MXDIM = 10000 )

      implicit real*8 (A-H,O-Z)

      character filen*20, textread*80

      dimension A(MXDIM)

      common /tapes/  IND

      common /cntrl/  alp, npu, npl, nnd
c
c---------------------------------------------------------------------c
c *** Problem definition :
c---------------------------------------------------------------------c
c
      IND = 51

      write(*,'(5(/),''    *** Program Loads 2D ***'',3(/))')
      filen = textread(' Enter the cp data file name : ')
c
c---------------------------------------------------------------------c
c *** Open input files :
c---------------------------------------------------------------------c
c
      write(*,'(/,'' >>> Opening data file ...'')')

      open(IND, file=filen, status='old', err=100)
c
c---------------------------------------------------------------------c
c *** Read problem size from top of data file :
c---------------------------------------------------------------------c
c
      write(*,'(/,'' >>> Getting problem size ...'')')

      rewind(IND)
      read(IND,*)
      read(IND,*) alp, npu, npl

	nnd = npu + npl

      write(*,'(''  >> No. of points along upper surface:'',I7)'),npu
      write(*,'(''  >> No. of points along lower surface:'',I7)'),npl
c
c---------------------------------------------------------------------c
c *** Create workspace pointers :
c---------------------------------------------------------------------c
c
      write(*,'(/,'' >>> Allocating workspace pointers ...'')')

      IEND = 1

      call ipoint(  IX,                 nnd, IEND )
      call ipoint(  IYP,                nnd, IEND )
      call ipoint(  ICP,                nnd, IEND )

      write(*,'(''  >> MEMORY ALLOCATED :     '', i10)') IEND
      write(*,'(''  >> MEMORY AVAILABLE :     '', i10)') MXDIM
      if ( IEND .GT. MXDIM ) then
         close(IND)
         stop 'Increase MXDIM !!!'
      endif
c
c---------------------------------------------------------------------c
c *** Read geometry data :
c---------------------------------------------------------------------c
c
      write(*,'(/,'' >>> Reading data file...'')')

      rewind(IND)
      call read_data( A(IX), A(IYP), A(ICP) )
      close(IND)
c
c---------------------------------------------------------------------c
c *** Compute the aerodynamic loads :
c---------------------------------------------------------------------c
c
      write(*,'(/,'' >>> Computing aerodynamic loads...'')')

      call compute_loads( A(IX), A(IYP), A(ICP) )
c
c---------------------------------------------------------------------c
c *** Clean up :
c---------------------------------------------------------------------c
c

c
c---------------------------------------------------------------------c
c *** Successful completion of program :
c---------------------------------------------------------------------c
c
  999 write (*,*) ' '
      stop ' OK!!'
c
c---------------------------------------------------------------------c
c *** Unsuccessful error messages :
c---------------------------------------------------------------------c
c
  100 write(*,'(/,'' Error opening data file.'')')
      stop
c
c---------------------------------------------------------------------c
c *** Format statements :
c---------------------------------------------------------------------c
c
  50  format( 3I5 )
c
      end
c*********************************************************************c
c*********************************************************************c




      character*80 function textread( prompt)
c*********************************************************************c
c
      character*(*) prompt
c
c---------------------------------------------------------------------c
c
      write(*,'(/,a,$)') prompt
      read(*,'(a)') textread
      return
c
      end
c*********************************************************************c
c*********************************************************************c




      integer function namlen( filen )
c*********************************************************************c
c
      character*20 filen
c
c---------------------------------------------------------------------c
c
      namlen = 0
      do i = 20,1,-1
          if ( filen(i:i) .ne. ' ' ) then
              namlen = i
              goto 101
          endif
      enddo
  101 return
c
      end
c*********************************************************************c
c*********************************************************************c





      subroutine ipoint( ipt, nsize, iend )
c*********************************************************************c
c
      ipt = iend
      iend = ipt + nsize
c
      return
c
      end
c*********************************************************************c
c*********************************************************************c





      subroutine read_data( X, YP, CP )
c*********************************************************************c
c
      implicit real*8 (A-H,O-Z)

      dimension X(nnd), YP(nnd), CP(nnd)

      common /tapes/  IND

      common /cntrl/  alp, npu, npl, nnd
c
c---------------------------------------------------------------------c
c *** Read the data file :
c---------------------------------------------------------------------c
c
c *** read the file header
c
      read(IND,*,err=100)
      read(IND,*,err=100) alp, npu, npl
      read(IND,*,err=100)

      if ( npu + npl .NE. nnd ) then
          close(IND)
          stop 'Unexpected error.'
      endif

c *** read the station coordinates, slope, and pressure coefficients
c
      do i = 1,nnd
          read(IND,*,err=100) X(i), YP(i), CP(i)
      enddo
c
c---------------------------------------------------------------------c
c *** Successful return :
c---------------------------------------------------------------------c
c
      return
c
c---------------------------------------------------------------------c
c *** Unsuccessful error messages :
c---------------------------------------------------------------------c
c
  100 write(*,'(/,'' Error reading data file.'')')
      close(IND)
      stop
c
      end
c*********************************************************************c
c*********************************************************************c









      subroutine compute_loads( X, YP, CP )
c*********************************************************************c
c
      implicit real*8 (A-H,O-Z)

      dimension X(nnd), YP(nnd), CP(nnd)

      common /cntrl/  alp, npu, npl, nnd
c
c---------------------------------------------------------------------c
c *** Initialize solution variables :
c---------------------------------------------------------------------c
c
      cx = 0.0d0
      cy = 0.0d0
      cm = 0.0d0

      pi = 4.0d0*DATAN(1.0d0)
c
c---------------------------------------------------------------------c
c *** Integrate forces along upper surface :
c---------------------------------------------------------------------c
c
      do i = 2,npu

c ******* define x & y-coords
c
          x1 =  X( i-1 )
          x2 =  X(  i  )
          y1 = YP( i-1 )
          y2 = YP(  i  )
	
c ******* define pressures
c
          p1 = CP( i-1 )
          p2 = CP(  i  )

c ******* accumulate forces
c
          cx = cx - 0.5d0*( p1*y1 + p2*y2 )*( x2 - x1 )
          cy = cy - 0.5d0*( p1 + p2 )*( x2 - x1 )

      enddo
c
c---------------------------------------------------------------------c
c *** Integrate forces along lower surface :
c---------------------------------------------------------------------c
c
      do i = 2,npl
	
c ******* define x & y-coords
c
          x1 =  X( npu + i-1 )
          x2 =  X(  npu + i  )
          y1 = YP( npu + i-1 )
          y2 = YP(  npu + i  )

c ******* define pressures
c
          p1 = CP( npu + i-1 )
          p2 = CP(  npu + i  )

c ******* accumulate forces
c
          cx = cx - 0.5d0*( p1*y1 + p2*y2 )*( x2 - x1 )
          cy = cy + 0.5d0*( p1 + p2 )*( x2 - x1 )

      enddo
c
c---------------------------------------------------------------------c
c *** Convert forces to aerodynamic loads :
c---------------------------------------------------------------------c
c
c *** convert angle of attack from degrees to radians
c
      alp = alp*pi/180.0d0

c *** rotate force vectors relative to free-stream
c
      cl = cy*DCOS(alp) - cx*DSIN(alp) 
      cd = cy*DSIN(alp) + cx*DCOS(alp) 
c
c---------------------------------------------------------------------c
c *** Output forces to screen :
c---------------------------------------------------------------------c
c
      write(*,'(''  >> Cx :     '', F10.5)') cx
      write(*,'(''  >> Cy :     '', F10.5)') cy

      write(*,*)
      write(*,'(''  >> Cl :     '', F10.5)') cl
      write(*,'(''  >> Cd :     '', F10.5)') cd
c      write(*,'(''  >> Cm :     '', F10.5)') cm
c
c---------------------------------------------------------------------c
c *** Successful return :
c---------------------------------------------------------------------c
c
      return
c
c---------------------------------------------------------------------c
c *** Unsuccessful error messages :
c---------------------------------------------------------------------c
c

c
      end
c*********************************************************************c
c*********************************************************************c