real*8 function vmarea ( betap, model )

C
C*******************************************************************************
C
C       Copyright (C) 1993, California Institute of Technology.  U.S.
C       Government Sponsorhip under NASA Contract NAS7-918 is
C       acknowledged.
C
C*******************************************************************************
C
C	Calculate the VMA Area of the Satellite
C
	double precision betap,vmatab(3,-90:90), betalo,betahi,alo,ahi,
     *			RESULT,x1,x2,xtemp, careas(2), breakpts(2)
	character model*4
	common /vma/vmatab, careas, breakpts

      common  /dragblock/  dragmodel, dragarea, atmos, atden
      character*10 dragmodel, atmos
      double precision dragarea, atden

        logical first
	data first/.true./

        IF (DRAGMODEL .EQ. 'CONSTANT') THEN
		VMARAEA = DRAGAREA
	        if (first) then
		   write(8,*) 'Using Constant Satellite Area.'
	           first = .false.
	        end if
		RETURN
	ELSE IF (DRAGMODEL .NE. 'VMA') THEN
		WRITE(8,*) 'INVALID DRAGMODEL = ', DRAGMODEL
		STOP '>>>>> ERROR EXIT.'
	END IF
        if (first) then
	   write(8,*) 'Using VMA Model for Satellite Area.'
           first = .false.
        end if

	if (model.eq.'DRAG') then
	  iarea = 2
	else if (model.eq.'SRP') then
          iarea = 3
        else
          write(8,*) 'Invalid VMA model = ',model
	  write(8,*) 'Valid models = DRAG, SRP'
	  stop '?? VMA Error Exit.'
	end if
C
C	interpolate between two successive values.  Find bounding values
C	in table.

	if ( model .eq. 'DRAG') then
	   if ( dabs(betap) .le. breakpts(1) ) then
	      result = careas(1)
	      goto 9999 
	   else if ( dabs(betap) .ge. breakpts(2) ) then
	      result = careas(2)
	      goto 9999 
	   end if
	end if

	jlow = -100

        do j=-90,89
	  x1 = vmatab(1,j)
          x2 = vmatab(1,j+1)
          if (x2.lt.x1) then
	    xtemp = x1
	    x1 = x2
            x2 = xtemp
	  end if

          if ( (x1.le.betap).and.(betap.le.x2)) then
	    jlow = j
	    jhi = j+1
	  end if

	end do

	if (jlow .lt. -90) then
	   write(8,*) 'Invalid Betap = ',betap
           stop 'Invalid Betap in VMA model.'
        end if
 
C
C	determine beta primes and areas for interpolation into table
C       
	betalo = vmatab(1,jlow)
	betahi = vmatab(1,jhi)	

	alow = vmatab(iarea, jlow)
        ahi  = vmatab(iarea, jhi)
C
C	linearly interpolate for the area
C
	RESULT = alow + (ahi-alow)*(betap-betalo)/(betahi-betalo)

C	WRITE(8,*) 'BETAP:',BETAP
C       WRITE(8,*) 'BETA BOUNDS:', BETAHI,BETALO
C	WRITE(8,*) 'BOUNDING AREAS:',ALOW,AHI
C	WRITE(8,*) 'INTERPOLATED AREA:',RESULT	
C

 9999   continue

C	write (8,*) betap,'  ', result

	VMAREA = RESULT
	return
	end