aboutsummaryrefslogtreecommitdiff
path: root/src/xgeom.f
diff options
context:
space:
mode:
Diffstat (limited to 'src/xgeom.f')
-rw-r--r--src/xgeom.f1794
1 files changed, 1794 insertions, 0 deletions
diff --git a/src/xgeom.f b/src/xgeom.f
new file mode 100644
index 0000000..c771632
--- /dev/null
+++ b/src/xgeom.f
@@ -0,0 +1,1794 @@
+C***********************************************************************
+C Module: xgeom.f
+C
+C Copyright (C) 2000 Mark Drela
+C
+C This program is free software; you can redistribute it and/or modify
+C it under the terms of the GNU General Public License as published by
+C the Free Software Foundation; either version 2 of the License, or
+C (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program; if not, write to the Free Software
+C Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+C***********************************************************************
+
+ SUBROUTINE LEFIND(SLE,X,XP,Y,YP,S,N)
+ DIMENSION X(*),XP(*),Y(*),YP(*),S(*)
+C------------------------------------------------------
+C Locates leading edge spline-parameter value SLE
+C
+C The defining condition is
+C
+C (X-XTE,Y-YTE) . (X',Y') = 0 at S = SLE
+C
+C i.e. the surface tangent is normal to the chord
+C line connecting X(SLE),Y(SLE) and the TE point.
+C------------------------------------------------------
+C
+C---- convergence tolerance
+ DSEPS = (S(N)-S(1)) * 1.0E-5
+C
+C---- set trailing edge point coordinates
+ XTE = 0.5*(X(1) + X(N))
+ YTE = 0.5*(Y(1) + Y(N))
+C
+C---- get first guess for SLE
+ DO 10 I=3, N-2
+ DXTE = X(I) - XTE
+ DYTE = Y(I) - YTE
+ DX = X(I+1) - X(I)
+ DY = Y(I+1) - Y(I)
+ DOTP = DXTE*DX + DYTE*DY
+ IF(DOTP .LT. 0.0) GO TO 11
+ 10 CONTINUE
+C
+ 11 SLE = S(I)
+C
+C---- check for sharp LE case
+ IF(S(I) .EQ. S(I-1)) THEN
+ccc WRITE(*,*) 'Sharp LE found at ',I,SLE
+ RETURN
+ ENDIF
+C
+C---- Newton iteration to get exact SLE value
+ DO 20 ITER=1, 50
+ XLE = SEVAL(SLE,X,XP,S,N)
+ YLE = SEVAL(SLE,Y,YP,S,N)
+ DXDS = DEVAL(SLE,X,XP,S,N)
+ DYDS = DEVAL(SLE,Y,YP,S,N)
+ DXDD = D2VAL(SLE,X,XP,S,N)
+ DYDD = D2VAL(SLE,Y,YP,S,N)
+C
+ XCHORD = XLE - XTE
+ YCHORD = YLE - YTE
+C
+C------ drive dot product between chord line and LE tangent to zero
+ RES = XCHORD*DXDS + YCHORD*DYDS
+ RESS = DXDS *DXDS + DYDS *DYDS
+ & + XCHORD*DXDD + YCHORD*DYDD
+C
+C------ Newton delta for SLE
+ DSLE = -RES/RESS
+C
+ DSLE = MAX( DSLE , -0.02*ABS(XCHORD+YCHORD) )
+ DSLE = MIN( DSLE , 0.02*ABS(XCHORD+YCHORD) )
+ SLE = SLE + DSLE
+ IF(ABS(DSLE) .LT. DSEPS) RETURN
+ 20 CONTINUE
+ WRITE(*,*) 'LEFIND: LE point not found. Continuing...'
+ SLE = S(I)
+ RETURN
+ END
+
+
+
+ SUBROUTINE XLFIND(SLE,X,XP,Y,YP,S,N)
+ DIMENSION X(*),XP(*),Y(*),YP(*),S(*)
+C------------------------------------------------------
+C Locates leftmost (minimum x) point location SLE
+C
+C The defining condition is
+C
+C X' = 0 at S = SLE
+C
+C i.e. the surface tangent is vertical
+C------------------------------------------------------
+C
+ DSLEN = S(N) - S(1)
+C
+C---- convergence tolerance
+ DSEPS = (S(N)-S(1)) * 1.0E-5
+C
+C---- get first guess for SLE
+ DO 10 I=3, N-2
+ DX = X(I+1) - X(I)
+ IF(DX .GT. 0.0) GO TO 11
+ 10 CONTINUE
+C
+ 11 SLE = S(I)
+C
+C---- check for sharp LE case
+ IF(S(I) .EQ. S(I-1)) THEN
+ccc WRITE(*,*) 'Sharp LE found at ',I,SLE
+ RETURN
+ ENDIF
+C
+C---- Newton iteration to get exact SLE value
+ DO 20 ITER=1, 50
+ DXDS = DEVAL(SLE,X,XP,S,N)
+ DXDD = D2VAL(SLE,X,XP,S,N)
+C
+C------ drive DXDS to zero
+ RES = DXDS
+ RESS = DXDD
+C
+C------ Newton delta for SLE
+ DSLE = -RES/RESS
+C
+ DSLE = MAX( DSLE , -0.01*ABS(DSLEN) )
+ DSLE = MIN( DSLE , 0.01*ABS(DSLEN) )
+ SLE = SLE + DSLE
+ IF(ABS(DSLE) .LT. DSEPS) RETURN
+ 20 CONTINUE
+ WRITE(*,*) 'XLFIND: Left point not found. Continuing...'
+ SLE = S(I)
+ RETURN
+ END ! XLFIND
+
+
+
+ SUBROUTINE NSFIND(SLE,X,XP,Y,YP,S,N)
+ REAL X(*),Y(*),S(*),XP(*),YP(*)
+C----------------------------------------------------------
+C Finds "nose" of airfoil where curvature is a maximum
+C----------------------------------------------------------
+C
+ PARAMETER (NMAX=500)
+ DIMENSION A(NMAX), B(NMAX), C(NMAX), CV(NMAX)
+C
+ IF(N.GT.NMAX) STOP 'NSFIND: Local array overflow. Increase NMAX.'
+C
+C---- set up curvature array
+ DO 3 I=1, N
+ CV(I) = CURV(S(I),X,XP,Y,YP,S,N)
+ 3 CONTINUE
+C
+C---- curvature smoothing length
+ SMOOL = 0.006*(S(N)-S(1))
+C
+C---- set up tri-diagonal system for smoothed curvatures
+ SMOOSQ = SMOOL**2
+ A(1) = 1.0
+ C(1) = 0.0
+ DO 4 I=2, N-1
+ DSM = S(I) - S(I-1)
+ DSP = S(I+1) - S(I)
+ DSO = 0.5*(S(I+1) - S(I-1))
+C
+ IF(DSM.EQ.0.0 .OR. DSP.EQ.0.0) THEN
+C------- leave curvature at corner point unchanged
+ B(I) = 0.0
+ A(I) = 1.0
+ C(I) = 0.0
+ ELSE
+ B(I) = SMOOSQ * ( - 1.0/DSM) / DSO
+ A(I) = SMOOSQ * ( 1.0/DSP + 1.0/DSM) / DSO + 1.0
+ C(I) = SMOOSQ * (-1.0/DSP ) / DSO
+ ENDIF
+ 4 CONTINUE
+ B(N) = 0.0
+ A(N) = 1.0
+C
+ CALL TRISOL(A,B,C,CV,N)
+C
+C---- find max curvature index
+ CVMAX = 0.
+ IVMAX = 0
+ DO 71 I=2, N-1
+ IF(ABS(CV(I)) .GT. CVMAX) THEN
+ CVMAX = ABS(CV(I))
+ IVMAX = I
+ ENDIF
+ 71 CONTINUE
+C
+C---- fit a parabola to the curvature at the three points near maximum
+ I = IVMAX
+C
+ IP = I+1
+ IM = I-1
+ IF(S(I) .EQ. S(IP)) IP = I+2
+ IF(S(I) .EQ. S(IM)) IM = I-2
+
+ DSM = S(I) - S(IM)
+ DSP = S(IP) - S(I)
+C
+ CVSM = (CV(I)-CV(IM))/DSM
+ CVSP = (CV(IP)-CV(I))/DSP
+C
+C---- 1st and 2nd derivatives at i=IVMAX
+ CVS = (CVSM*DSP + CVSP*DSM)/(DSP+DSM)
+ CVSS = 2.0*(CVSP-CVSM)/(DSP+DSM)
+C
+C---- set location of arc length at maximum of parabola
+ DS = -CVS/CVSS
+ SLE = S(I) + DS
+C
+ RETURN
+ END
+
+
+ SUBROUTINE SOPPS(SOPP, SI, X,XP,Y,YP,S,N, SLE)
+ DIMENSION X(*),XP(*),Y(*),YP(*),S(*)
+C--------------------------------------------------
+C Calculates arc length SOPP of point
+C which is opposite of point SI, on the
+C other side of the airfoil baseline
+C--------------------------------------------------
+C
+C---- reference length for testing convergence
+ SLEN = S(N) - S(1)
+C
+C---- set chordline vector
+ XLE = SEVAL(SLE,X,XP,S,N)
+ YLE = SEVAL(SLE,Y,YP,S,N)
+ XTE = 0.5*(X(1)+X(N))
+ YTE = 0.5*(Y(1)+Y(N))
+ CHORD = SQRT((XTE-XLE)**2 + (YTE-YLE)**2)
+ DXC = (XTE-XLE) / CHORD
+ DYC = (YTE-YLE) / CHORD
+C
+ IF(SI.LT.SLE) THEN
+ IN = 1
+ INOPP = N
+ ELSE
+ IN = N
+ INOPP = 1
+ ENDIF
+ SFRAC = (SI-SLE)/(S(IN)-SLE)
+ SOPP = SLE + SFRAC*(S(INOPP)-SLE)
+C
+ IF(ABS(SFRAC) .LE. 1.0E-5) THEN
+ SOPP = SLE
+ RETURN
+ ENDIF
+C
+C---- XBAR = x coordinate in chord-line axes
+ XI = SEVAL(SI , X,XP,S,N)
+ YI = SEVAL(SI , Y,YP,S,N)
+ XLE = SEVAL(SLE, X,XP,S,N)
+ YLE = SEVAL(SLE, Y,YP,S,N)
+ XBAR = (XI-XLE)*DXC + (YI-YLE)*DYC
+C
+C---- converge on exact opposite point with same XBAR value
+ DO 300 ITER=1, 12
+ XOPP = SEVAL(SOPP,X,XP,S,N)
+ YOPP = SEVAL(SOPP,Y,YP,S,N)
+ XOPPD = DEVAL(SOPP,X,XP,S,N)
+ YOPPD = DEVAL(SOPP,Y,YP,S,N)
+C
+ RES = (XOPP -XLE)*DXC + (YOPP -YLE)*DYC - XBAR
+ RESD = XOPPD *DXC + YOPPD *DYC
+C
+ IF(ABS(RES)/SLEN .LT. 1.0E-5) GO TO 305
+ IF(RESD .EQ. 0.0) GO TO 303
+C
+ DSOPP = -RES/RESD
+ SOPP = SOPP + DSOPP
+C
+ IF(ABS(DSOPP)/SLEN .LT. 1.0E-5) GO TO 305
+ 300 CONTINUE
+ 303 WRITE(*,*)
+ & 'SOPPS: Opposite-point location failed. Continuing...'
+ SOPP = SLE + SFRAC*(S(INOPP)-SLE)
+C
+ 305 CONTINUE
+ RETURN
+ END ! SOPPS
+
+
+
+ SUBROUTINE NORM(X,XP,Y,YP,S,N)
+ DIMENSION X(*),XP(*),Y(*),YP(*),S(*)
+C-----------------------------------------------
+C Scales coordinates to get unit chord
+C-----------------------------------------------
+C
+ CALL SCALC(X,Y,S,N)
+ CALL SEGSPL(X,XP,S,N)
+ CALL SEGSPL(Y,YP,S,N)
+C
+ CALL LEFIND(SLE,X,XP,Y,YP,S,N)
+C
+ XMAX = 0.5*(X(1) + X(N))
+ XMIN = SEVAL(SLE,X,XP,S,N)
+ YMIN = SEVAL(SLE,Y,YP,S,N)
+C
+ FUDGE = 1.0/(XMAX-XMIN)
+ DO 40 I=1, N
+ X(I) = (X(I)-XMIN)*FUDGE
+ Y(I) = (Y(I)-YMIN)*FUDGE
+ S(I) = S(I)*FUDGE
+ 40 CONTINUE
+C
+ RETURN
+ END
+
+
+ SUBROUTINE GEOPAR(X,XP,Y,YP,S,N, T,
+ & SLE,CHORD,AREA,RADLE,ANGTE,
+ & EI11A,EI22A,APX1A,APX2A,
+ & EI11T,EI22T,APX1T,APX2T,
+ & THICK,CAMBR)
+ DIMENSION X(*), XP(*), Y(*), YP(*), S(*), T(*)
+C
+ PARAMETER (IBX=600)
+ DIMENSION
+ & XCAM(2*IBX), YCAM(2*IBX), YCAMP(2*IBX),
+ & XTHK(2*IBX), YTHK(2*IBX), YTHKP(2*IBX)
+C------------------------------------------------------
+C Sets geometric parameters for airfoil shape
+C------------------------------------------------------
+ CALL LEFIND(SLE,X,XP,Y,YP,S,N)
+C
+ XLE = SEVAL(SLE,X,XP,S,N)
+ YLE = SEVAL(SLE,Y,YP,S,N)
+ XTE = 0.5*(X(1)+X(N))
+ YTE = 0.5*(Y(1)+Y(N))
+C
+ CHSQ = (XTE-XLE)**2 + (YTE-YLE)**2
+ CHORD = SQRT(CHSQ)
+C
+ CURVLE = CURV(SLE,X,XP,Y,YP,S,N)
+C
+ RADLE = 0.0
+ IF(ABS(CURVLE) .GT. 0.001*(S(N)-S(1))) RADLE = 1.0 / CURVLE
+C
+ ANG1 = ATAN2( -YP(1) , -XP(1) )
+ ANG2 = ATANC( YP(N) , XP(N) , ANG1 )
+ ANGTE = ANG2 - ANG1
+C
+
+ DO I=1, N
+ T(I) = 1.0
+ ENDDO
+C
+ CALL AECALC(N,X,Y,T, 1,
+ & AREA,XCENA,YCENA,EI11A,EI22A,APX1A,APX2A)
+C
+ CALL AECALC(N,X,Y,T, 2,
+ & SLEN,XCENT,YCENT,EI11T,EI22T,APX1T,APX2T)
+C
+C--- Old, approximate thickness,camber routine (on discrete points only)
+ CALL TCCALC(X,XP,Y,YP,S,N, THICK,XTHICK, CAMBR,XCAMBR )
+C
+C--- More accurate thickness and camber estimates
+cc CALL GETCAM(XCAM,YCAM,NCAM,XTHK,YTHK,NTHK,
+cc & X,XP,Y,YP,S,N )
+cc CALL GETMAX(XCAM,YCAM,YCAMP,NCAM,XCAMBR,CAMBR)
+cc CALL GETMAX(XTHK,YTHK,YTHKP,NTHK,XTHICK,THICK)
+cc THICK = 2.0*THICK
+C
+ WRITE(*,1000) THICK,XTHICK,CAMBR,XCAMBR
+ 1000 FORMAT( ' Max thickness = ',F12.6,' at x = ',F7.3,
+ & /' Max camber = ',F12.6,' at x = ',F7.3)
+
+
+C
+ RETURN
+ END ! GEOPAR
+
+
+ SUBROUTINE AECALC(N,X,Y,T, ITYPE,
+ & AREA,XCEN,YCEN,EI11,EI22,APX1,APX2)
+ DIMENSION X(*),Y(*),T(*)
+C---------------------------------------------------------------
+C Calculates geometric properties of shape X,Y
+C
+C Input:
+C N number of points
+C X(.) shape coordinate point arrays
+C Y(.)
+C T(.) skin-thickness array, used only if ITYPE = 2
+C ITYPE = 1 ... integration is over whole area dx dy
+C = 2 ... integration is over skin area t ds
+C
+C Output:
+C XCEN,YCEN centroid location
+C EI11,EI22 principal moments of inertia
+C APX1,APX2 principal-axis angles
+C---------------------------------------------------------------
+ DATA PI / 3.141592653589793238 /
+C
+ SINT = 0.0
+ AINT = 0.0
+ XINT = 0.0
+ YINT = 0.0
+ XXINT = 0.0
+ XYINT = 0.0
+ YYINT = 0.0
+C
+ DO 10 IO = 1, N
+ IF(IO.EQ.N) THEN
+ IP = 1
+ ELSE
+ IP = IO + 1
+ ENDIF
+C
+ DX = X(IO) - X(IP)
+ DY = Y(IO) - Y(IP)
+ XA = (X(IO) + X(IP))*0.50
+ YA = (Y(IO) + Y(IP))*0.50
+ TA = (T(IO) + T(IP))*0.50
+C
+ DS = SQRT(DX*DX + DY*DY)
+ SINT = SINT + DS
+
+ IF(ITYPE.EQ.1) THEN
+C-------- integrate over airfoil cross-section
+ DA = YA*DX
+ AINT = AINT + DA
+ XINT = XINT + XA *DA
+ YINT = YINT + YA *DA/2.0
+ XXINT = XXINT + XA*XA*DA
+ XYINT = XYINT + XA*YA*DA/2.0
+ YYINT = YYINT + YA*YA*DA/3.0
+ ELSE
+C-------- integrate over skin thickness
+ DA = TA*DS
+ AINT = AINT + DA
+ XINT = XINT + XA *DA
+ YINT = YINT + YA *DA
+ XXINT = XXINT + XA*XA*DA
+ XYINT = XYINT + XA*YA*DA
+ YYINT = YYINT + YA*YA*DA
+ ENDIF
+C
+ 10 CONTINUE
+C
+ AREA = AINT
+C
+ IF(AINT .EQ. 0.0) THEN
+ XCEN = 0.0
+ YCEN = 0.0
+ EI11 = 0.0
+ EI22 = 0.0
+ APX1 = 0.0
+ APX2 = ATAN2(1.0,0.0)
+ RETURN
+ ENDIF
+C
+C
+C---- calculate centroid location
+ XCEN = XINT/AINT
+ YCEN = YINT/AINT
+C
+C---- calculate inertias
+ EIXX = YYINT - YCEN*YCEN*AINT
+ EIXY = XYINT - XCEN*YCEN*AINT
+ EIYY = XXINT - XCEN*XCEN*AINT
+C
+C---- set principal-axis inertias, EI11 is closest to "up-down" bending inertia
+ EISQ = 0.25*(EIXX - EIYY)**2 + EIXY**2
+ SGN = SIGN( 1.0 , EIYY-EIXX )
+ EI11 = 0.5*(EIXX + EIYY) - SGN*SQRT(EISQ)
+ EI22 = 0.5*(EIXX + EIYY) + SGN*SQRT(EISQ)
+C
+ IF(EI11.EQ.0.0 .OR. EI22.EQ.0.0) THEN
+C----- vanishing section stiffness
+ APX1 = 0.0
+ APX2 = ATAN2(1.0,0.0)
+C
+ ELSEIF(EISQ/(EI11*EI22) .LT. (0.001*SINT)**4) THEN
+C----- rotationally-invariant section (circle, square, etc.)
+ APX1 = 0.0
+ APX2 = ATAN2(1.0,0.0)
+C
+ ELSE
+C----- normal airfoil section
+ C1 = EIXY
+ S1 = EIXX-EI11
+C
+ C2 = EIXY
+ S2 = EIXX-EI22
+C
+ IF(ABS(S1).GT.ABS(S2)) THEN
+ APX1 = ATAN2(S1,C1)
+ APX2 = APX1 + 0.5*PI
+ ELSE
+ APX2 = ATAN2(S2,C2)
+ APX1 = APX2 - 0.5*PI
+ ENDIF
+
+ IF(APX1.LT.-0.5*PI) APX1 = APX1 + PI
+ IF(APX1.GT.+0.5*PI) APX1 = APX1 - PI
+ IF(APX2.LT.-0.5*PI) APX2 = APX2 + PI
+ IF(APX2.GT.+0.5*PI) APX2 = APX2 - PI
+C
+ ENDIF
+C
+ RETURN
+ END ! AECALC
+
+
+
+ SUBROUTINE TCCALC(X,XP,Y,YP,S,N,
+ & THICK,XTHICK, CAMBR,XCAMBR )
+ DIMENSION X(*),XP(*),Y(*),YP(*),S(*)
+C---------------------------------------------------------------
+C Calculates max thickness and camber at airfoil points
+C
+C Note: this routine does not find the maximum camber or
+C thickness exactly as it only looks at discrete points
+C
+C Input:
+C N number of points
+C X(.) shape coordinate point arrays
+C Y(.)
+C
+C Output:
+C THICK max thickness
+C CAMBR max camber
+C---------------------------------------------------------------
+ CALL LEFIND(SLE,X,XP,Y,YP,S,N)
+ XLE = SEVAL(SLE,X,XP,S,N)
+ YLE = SEVAL(SLE,Y,YP,S,N)
+ XTE = 0.5*(X(1)+X(N))
+ YTE = 0.5*(Y(1)+Y(N))
+ CHORD = SQRT((XTE-XLE)**2 + (YTE-YLE)**2)
+C
+C---- set unit chord-line vector
+ DXC = (XTE-XLE) / CHORD
+ DYC = (YTE-YLE) / CHORD
+C
+ THICK = 0.
+ XTHICK = 0.
+ CAMBR = 0.
+ XCAMBR = 0.
+C
+C---- go over each point, finding the y-thickness and camber
+ DO 30 I=1, N
+ XBAR = (X(I)-XLE)*DXC + (Y(I)-YLE)*DYC
+ YBAR = (Y(I)-YLE)*DXC - (X(I)-XLE)*DYC
+C
+C------ set point on the opposite side with the same chord x value
+ CALL SOPPS(SOPP, S(I), X,XP,Y,YP,S,N, SLE)
+ XOPP = SEVAL(SOPP,X,XP,S,N)
+ YOPP = SEVAL(SOPP,Y,YP,S,N)
+C
+ YBAROP = (YOPP-YLE)*DXC - (XOPP-XLE)*DYC
+C
+ YC = 0.5*(YBAR+YBAROP)
+ YT = ABS(YBAR-YBAROP)
+C
+ IF(ABS(YC) .GT. ABS(CAMBR)) THEN
+ CAMBR = YC
+ XCAMBR = XOPP
+ ENDIF
+ IF(ABS(YT) .GT. ABS(THICK)) THEN
+ THICK = YT
+ XTHICK = XOPP
+ ENDIF
+ 30 CONTINUE
+C
+ RETURN
+ END ! TCCALC
+
+
+
+ SUBROUTINE YSYM(X,XP,Y,YP,S,NX,N,ISIDE, XNEW,YNEW)
+C---------------------------------------------------------
+C Makes passed-in airfoil symmetric about chord line.
+C---------------------------------------------------------
+
+ DIMENSION X(NX),XP(NX),Y(NX),YP(NX),S(NX)
+ DIMENSION XNEW(NX), YNEW(NX)
+C
+ SREF = S(N) - S(1)
+C
+ CALL LEFIND(SLE,X,XP,Y,YP,S,N)
+ XLE = SEVAL(SLE,X,XP,S,N)
+ YLE = SEVAL(SLE,Y,YP,S,N)
+ XTE = 0.5*(X(1)+X(N))
+ YTE = 0.5*(Y(1)+Y(N))
+ CHSQ = (XTE-XLE)**2 + (YTE-YLE)**2
+C
+C---- set unit chord-line vector
+ DXC = (XTE-XLE) / SQRT(CHSQ)
+ DYC = (YTE-YLE) / SQRT(CHSQ)
+C
+C---- find index of node ILE which is just before leading edge point
+ DO 5 I=2, N
+ DS = S(I) - S(I-1)
+ IF(S(I)-SLE .GE. -0.01*DS) GO TO 6
+ 5 CONTINUE
+ 6 CONTINUE
+ ILE = I-1
+C
+ DS = S(ILE+1) - S(ILE)
+ IF(SLE-S(ILE-1) .LT. 0.1*DS) THEN
+C------ point is just before LE, we will move it ahead to LE
+ ILE1 = ILE - 1
+ ILE2 = ILE + 1
+ ELSE IF(S(ILE+1)-SLE .LT. 0.1*DS) THEN
+C------ point is just after LE, we will move it back to LE
+ ILE1 = ILE
+ ILE2 = ILE + 2
+ ELSE
+C------ no point is near LE ... we will add new point
+ ILE1 = ILE
+ ILE2 = ILE + 1
+ ENDIF
+C
+C---- set index limits of side which will set symmetric geometry
+ IF(ISIDE.EQ.1) THEN
+ IG1 = 1
+ IG2 = ILE1
+ IGDIR = +1
+ ELSE
+ IG1 = N
+ IG2 = ILE2
+ IGDIR = -1
+ ENDIF
+C
+C---- set new number of points, including LE point
+ NNEW = 2*(IABS(IG2-IG1) + 1) + 1
+ IF(NNEW.GT.NX) STOP 'YSYM: Array overflow on passed arrays.'
+C
+C---- set symmetric geometry
+ DO 10 I=IG1, IG2, IGDIR
+C
+C------ coordinates in chord-line axes
+ XBAR = (X(I)-XLE)*DXC + (Y(I)-YLE)*DYC
+ YBAR = (Y(I)-YLE)*DXC - (X(I)-XLE)*DYC
+C
+ I1 = 1 + (I - IG1)*IGDIR
+ I2 = NNEW - (I - IG1)*IGDIR
+C
+ XNEW(I1) = XLE + XBAR*DXC - YBAR*DYC
+ XNEW(I2) = XLE + XBAR*DXC + YBAR*DYC
+C
+ YNEW(I1) = YLE + YBAR*DXC + XBAR*DYC
+ YNEW(I2) = YLE - YBAR*DXC + XBAR*DYC
+ 10 CONTINUE
+C
+C---- set new LE point
+ XNEW(NNEW/2+1) = XLE
+ YNEW(NNEW/2+1) = YLE
+C
+C---- set geometry for returning
+ N = NNEW
+ DO 20 IG = 1, N
+ IF(IGDIR.EQ.+1) THEN
+ I = IG
+ ELSE
+ I = N - IG + 1
+ ENDIF
+ X(I) = XNEW(IG)
+ Y(I) = YNEW(IG)
+ 20 CONTINUE
+C
+ CALL SCALC(X,Y,S,N)
+ CALL SEGSPL(X,XP,S,N)
+ CALL SEGSPL(Y,YP,S,N)
+C
+ RETURN
+ END ! YSYM
+
+
+
+ SUBROUTINE LERSCL(X,XP,Y,YP,S,N, DOC,RFAC, XNEW,YNEW)
+C---------------------------------------------------------
+C Adjusts airfoil to scale LE radius by factor RFAC.
+C Blending of new shape is done with decay length DOC.
+C---------------------------------------------------------
+ DIMENSION X(*),XP(*),Y(*),YP(*),S(*)
+ DIMENSION XNEW(*), YNEW(*)
+C
+ CALL LEFIND(SLE,X,XP,Y,YP,S,N)
+ XLE = SEVAL(SLE,X,XP,S,N)
+ YLE = SEVAL(SLE,Y,YP,S,N)
+ XTE = 0.5*(X(1)+X(N))
+ YTE = 0.5*(Y(1)+Y(N))
+ CHORD = SQRT((XTE-XLE)**2 + (YTE-YLE)**2)
+C
+C---- set unit chord-line vector
+ DXC = (XTE-XLE) / CHORD
+ DYC = (YTE-YLE) / CHORD
+C
+ SRFAC = SQRT(ABS(RFAC))
+C
+C---- go over each point, changing the y-thickness appropriately
+ DO 30 I=1, N
+ XBAR = (X(I)-XLE)*DXC + (Y(I)-YLE)*DYC
+ YBAR = (Y(I)-YLE)*DXC - (X(I)-XLE)*DYC
+C
+C------ set point on the opposite side with the same chord x value
+ CALL SOPPS(SOPP, S(I), X,XP,Y,YP,S,N, SLE)
+ XOPP = SEVAL(SOPP,X,XP,S,N)
+ YOPP = SEVAL(SOPP,Y,YP,S,N)
+C
+ YBAROP = (YOPP-YLE)*DXC - (XOPP-XLE)*DYC
+C
+C------ thickness factor tails off exponentially towards trailing edge
+ XOC = XBAR/CHORD
+ ARG = MIN( XOC/DOC , 15.0 )
+ TFAC = 1.0 - (1.0-SRFAC)*EXP(-ARG)
+C
+C------ set new chord x,y coordinates by changing thickness locally
+ YBARCT = 0.5*(YBAR+YBAROP) + TFAC*0.5*(YBAR-YBAROP)
+C
+ XNEW(I) = XLE + XBAR *DXC - YBARCT*DYC
+ YNEW(I) = YLE + YBARCT*DXC + XBAR *DYC
+ 30 CONTINUE
+C
+ RETURN
+ END
+
+
+
+ SUBROUTINE SSS(SS,S1,S2,DEL,XBF,YBF,X,XP,Y,YP,S,N,ISIDE)
+ DIMENSION X(*),XP(*),Y(*),YP(*),S(*)
+C----------------------------------------------------------------
+C Returns arc length points S1,S2 at flap surface break
+C locations. S1 is on fixed airfoil part, S2 is on flap.
+C The points are defined according to two cases:
+C
+C
+C If DEL > 0: Surface will be eliminated in S1 < s < S2
+C
+C Returns the arc length values S1,S2 of the endpoints
+C of the airfoil surface segment which "disappears" as a
+C result of the flap deflection. The line segments between
+C these enpoints and the flap hinge point (XBF,YBF) have
+C an included angle of DEL. DEL is therefore the flap
+C deflection which will join up the points at S1,S2.
+C SS is an approximate arc length value near S1 and S2.
+C It is used as an initial guess for the Newton loop
+C for S1 and S2.
+C
+C
+C If DEL = 0: Surface will be created at s = S1 = S2
+C
+C If DEL=0, then S1,S2 will cooincide, and will be located
+C on the airfoil surface where the segment joining the
+C point at S1,S2 and the hinge point is perpendicular to
+C the airfoil surface. This will be the point where the
+C airfoil surface must be broken to permit a gap to open
+C as a result of the flap deflection.
+C----------------------------------------------------------------
+C
+C---- convergence epsilon
+ DATA EPS / 1.0E-5 /
+C
+ STOT = ABS( S(N) - S(1) )
+C
+ SIND = SIN(0.5*ABS(DEL))
+C
+ SSGN = 1.0
+ IF(ISIDE.EQ.1) SSGN = -1.0
+C
+C---- initial guesses for S1, S2
+ RSQ = (SEVAL(SS,X,XP,S,N)-XBF)**2 + (SEVAL(SS,Y,YP,S,N)-YBF)**2
+ S1 = SS - (SIND*SQRT(RSQ) + EPS*STOT)*SSGN
+ S2 = SS + (SIND*SQRT(RSQ) + EPS*STOT)*SSGN
+C
+C---- Newton iteration loop
+ DO 10 ITER=1, 10
+ X1 = SEVAL(S1,X,XP,S,N)
+ X1P = DEVAL(S1,X,XP,S,N)
+ Y1 = SEVAL(S1,Y,YP,S,N)
+ Y1P = DEVAL(S1,Y,YP,S,N)
+C
+ X2 = SEVAL(S2,X,XP,S,N)
+ X2P = DEVAL(S2,X,XP,S,N)
+ Y2 = SEVAL(S2,Y,YP,S,N)
+ Y2P = DEVAL(S2,Y,YP,S,N)
+C
+ R1SQ = (X1-XBF)**2 + (Y1-YBF)**2
+ R2SQ = (X2-XBF)**2 + (Y2-YBF)**2
+ R1 = SQRT(R1SQ)
+ R2 = SQRT(R2SQ)
+C
+ RRSQ = (X1-X2)**2 + (Y1-Y2)**2
+ RR = SQRT(RRSQ)
+C
+ IF(R1.LE.EPS*STOT .OR. R2.LE.EPS*STOT) THEN
+ S1 = SS
+ S2 = SS
+ RETURN
+ ENDIF
+C
+ R1_S1 = (X1P*(X1-XBF) + Y1P*(Y1-YBF))/R1
+ R2_S2 = (X2P*(X2-XBF) + Y2P*(Y2-YBF))/R2
+C
+ IF(SIND.GT.0.01) THEN
+C
+ IF(RR.EQ.0.0) RETURN
+C
+ RR_S1 = (X1P*(X1-X2) + Y1P*(Y1-Y2))/RR
+ RR_S2 = -(X2P*(X1-X2) + Y2P*(Y1-Y2))/RR
+C
+C------- Residual 1: set included angle via dot product
+ RS1 = ((XBF-X1)*(X2-X1) + (YBF-Y1)*(Y2-Y1))/RR - SIND*R1
+ A11 = ((XBF-X1)*( -X1P) + (YBF-Y1)*( -Y1P))/RR
+ & + (( -X1P)*(X2-X1) + ( -Y1P)*(Y2-Y1))/RR
+ & - ((XBF-X1)*(X2-X1) + (YBF-Y1)*(Y2-Y1))*RR_S1/RRSQ
+ & - SIND*R1_S1
+ A12 = ((XBF-X1)*(X2P ) + (YBF-Y1)*(Y2P ))/RR
+ & - ((XBF-X1)*(X2-X1) + (YBF-Y1)*(Y2-Y1))*RR_S2/RRSQ
+C
+C------- Residual 2: set equal length segments
+ RS2 = R1 - R2
+ A21 = R1_S1
+ A22 = - R2_S2
+C
+ ELSE
+C
+C------- Residual 1: set included angle via small angle approximation
+ RS1 = (R1+R2)*SIND + (S1 - S2)*SSGN
+ A11 = R1_S1 *SIND + SSGN
+ A12 = R2_S2 *SIND - SSGN
+C
+C------- Residual 2: set vector sum of line segments beteen the
+C- endpoints and flap hinge to be perpendicular to airfoil surface.
+ X1PP = D2VAL(S1,X,XP,S,N)
+ Y1PP = D2VAL(S1,Y,YP,S,N)
+ X2PP = D2VAL(S2,X,XP,S,N)
+ Y2PP = D2VAL(S2,Y,YP,S,N)
+C
+ XTOT = X1+X2 - 2.0*XBF
+ YTOT = Y1+Y2 - 2.0*YBF
+C
+ RS2 = XTOT*(X1P+X2P) + YTOT*(Y1P+Y2P)
+ A21 = X1P*(X1P+X2P) + Y1P*(Y1P+Y2P) + XTOT*X1PP + YTOT*Y1PP
+ A22 = X2P*(X1P+X2P) + Y2P*(Y1P+Y2P) + XTOT*X2PP + YTOT*Y2PP
+C
+ ENDIF
+C
+ DET = A11*A22 - A12*A21
+ DS1 = -(RS1*A22 - A12*RS2) / DET
+ DS2 = -(A11*RS2 - RS1*A21) / DET
+C
+ DS1 = MIN( DS1 , 0.01*STOT )
+ DS1 = MAX( DS1 , -.01*STOT )
+ DS2 = MIN( DS2 , 0.01*STOT )
+ DS2 = MAX( DS2 , -.01*STOT )
+C
+ S1 = S1 + DS1
+ S2 = S2 + DS2
+ IF(ABS(DS1)+ABS(DS2) .LT. EPS*STOT ) GO TO 11
+ 10 CONTINUE
+ WRITE(*,*) 'SSS: failed to converge subtending angle points'
+ S1 = SS
+ S2 = SS
+C
+ 11 CONTINUE
+C
+C---- make sure points are identical if included angle is zero.
+ IF(DEL.EQ.0.0) THEN
+ S1 = 0.5*(S1+S2)
+ S2 = S1
+ ENDIF
+C
+ RETURN
+ END
+
+
+ SUBROUTINE CLIS(X,XP,Y,YP,S,N)
+ DIMENSION X(*), XP(*), Y(*), YP(*), S(*)
+C-------------------------------------------------------------------
+C Displays curvatures at panel nodes.
+C-------------------------------------------------------------------
+ PI = 4.0*ATAN(1.0)
+C
+ CMAX = 0.0
+ IMAX = 1
+C
+C---- go over each point, calculating curvature
+ WRITE(*,1050)
+ DO 30 I=1, N
+ IF(I.EQ.1) THEN
+ ARAD = ATAN2(-YP(I),-XP(I))
+ ELSE
+ ARAD = ATANC(-YP(I),-XP(I),ARAD)
+ ENDIF
+ ADEG = ARAD * 180.0/PI
+ CV = CURV(S(I),X,XP,Y,YP,S,N)
+ WRITE(*,1100) I, X(I), Y(I), S(I), ADEG, CV
+ IF(ABS(CV) .GT. ABS(CMAX)) THEN
+ CMAX = CV
+ IMAX = I
+ ENDIF
+ 30 CONTINUE
+C
+ WRITE(*,1200) CMAX, IMAX, X(IMAX), Y(IMAX), S(IMAX)
+C
+ RETURN
+C
+ 1050 FORMAT(
+ & /' i x y s theta curv')
+CCC 120 0.12134 -0.10234 -0.30234 180.024 2025.322
+ 1100 FORMAT(1X,I3, 3F10.5, F11.3, F12.3)
+ 1200 FORMAT(/' Maximum curvature =', F14.3,
+ & ' at i,x,y,s = ', I3, 3F9.4 )
+ END ! CLIS
+
+
+
+
+ SUBROUTINE PLTCRV(SBLE,XB,XBP,YB,YBP,SB,NB,CV)
+C
+C---- Plot the curvature on the blade
+C
+ DIMENSION XB(NB),XBP(NB),YB(NB),YBP(NB),SB(NB),CV(NB)
+ CHARACTER ANS*1, ANSARG*128
+ LOGICAL LCVEXP, ERROR
+C
+ DATA LMASK0, LMASK1, LMASK2, LMASK3 / -1, -32640, -30584, -21846 /
+C
+ CH = 0.01
+ LCVEXP = .FALSE.
+ CVEXP = 1.0/3.0
+C
+ 10 SBTOT = 0.5*(SB(NB)-SB(1))
+ XTE = 0.5*(XB(NB)+XB(1))
+ YTE = 0.5*(YB(NB)+YB(1))
+ XLE = SEVAL(SBLE,XB,XBP,SB,NB)
+ YLE = SEVAL(SBLE,YB,YBP,SB,NB)
+ CVLE = CURV(SBLE,XB,XBP,YB,YBP,SB,NB) * SBTOT
+ IF(LCVEXP) CVLE = CVLE**CVEXP
+C
+ CVMAX = CVLE
+ SVMAX = SBLE
+ XVMAX = XLE
+ CVMIN = CVLE
+ SVMIN = SBLE
+ XVMIN = XLE
+ DO I=1, NB
+C---- set up curvature array
+ CV(I) = CURV(SB(I),XB,XBP,YB,YBP,SB,NB) * SBTOT
+ IF(LCVEXP) THEN
+ IF(CV(I).GT.0.0) THEN
+ CV(I) = CV(I)**CVEXP
+ ELSEIF(CV(I).EQ.0.0) THEN
+ CV(I) = 0.0
+ ELSEIF(CV(I).LT.0.0) THEN
+ CVSGN = SIGN(1.0,CV(I))
+ CV(I) = CVSGN*(ABS(CV(I))**CVEXP)
+ ENDIF
+ ENDIF
+ IF(CV(I).GT.CVMAX) THEN
+ CVMAX = CV(I)
+ SVMAX = SB(I)
+ XVMAX = XB(I)
+ ENDIF
+ IF(CV(I).LT.CVMIN) THEN
+ CVMIN = CV(I)
+ SVMIN = SB(I)
+ XVMIN = XB(I)
+ ENDIF
+ IF(SB(I).LE.SBLE) ILE = I
+ END DO
+C
+cc CALL SCALIT(1,CVMAX-CVMIN,0.0,CWT)
+ CVMX = CVMAX
+ CVMN = CVMIN
+ CALL AXISADJ(CVMN,CVMX,CVSPAN,CVDEL,NCVTICS)
+ CWT = 1.0/CVSPAN
+ XMX = XTE
+ XMN = XLE
+ CALL AXISADJ(XMN,XMX,XSPAN,XDEL,NXTICS)
+C--- Correct min/max for points just slightly off from a major division in X
+ IF(XLE-XMN.GT.0.95*XDEL) XMN = XMN + XDEL
+ IF(XMX-XTE.GT.0.95*XDEL) XMX = XMX - XDEL
+ XSPAN = XMX-XMN
+ XWT = 1.0/XSPAN
+C
+ PAR = 0.75
+ XLEN = 0.8
+ YLEN = PAR*XLEN
+ XMIN = XMN
+ XMAX = XMX
+ XDEL = XDEL
+ NXG = (XMAX-XMIN)/XDEL
+ XSF = XLEN/(XMAX-XMIN)
+ XOF = XMN
+ YMIN = CVMN
+ YMAX = CVMX
+ YDEL = CVDEL
+ NYG = (YMAX-YMIN)/YDEL
+ YSF = YLEN/(YMAX-YMIN)
+ YOF = 0.0
+C
+ CALL PLTINI
+ CALL PLOT(0.14,0.1+YLEN*(-YMIN/(YMAX-YMIN)),-3)
+C
+C--- X axis (x/c)
+ CALL NEWPEN(2)
+ XLN = XLEN
+ IF(XMIN.EQ.0.0) XLN = -XLN
+ CALL XAXIS(0.0,0.0,XLN,XSF*XDEL,XMIN,XDEL,CH,1)
+ XC = XSF*3.5*XDEL -0.5*1.2*CH
+ YC = -3.5*1.2*CH
+ CALL PLCHAR(XC,YC,1.2*CH,'X',0.0,1)
+C
+C--- Y axis (curvature)
+ CALL YAXIS(0.0,YSF*YMIN,YLEN,YSF*YDEL,YMIN,YDEL,CH,1)
+ XC = -4.5*1.2*CH
+ YC = YSF*(YMAX-0.5*YDEL) - 0.5*1.2*CH
+ IF(LCVEXP) THEN
+ CALL PLCHAR(XC-4.5*1.2*CH,YC,1.2*CH,'CV^',0.0,3)
+ CALL PLNUMB(XC-1.5*1.2*CH,YC+0.5*CH,CH,CVEXP,0.0,2)
+ ELSE
+ CALL PLCHAR(XC,YC,1.2*CH,'CV',0.0,2)
+ ENDIF
+C
+ CALL PLGRID(0.0,YSF*YMIN, NXG,XSF*XDEL, NYG,YSF*YDEL, LMASK2)
+ XC = 0.0
+ YC = YSF*YMAX + 1.0*1.2*CH
+ IF(LCVEXP) THEN
+ CALL PLCHAR(XC,YC,1.2*CH,'Curvature^n vs X',0.0,16)
+ ELSE
+ CALL PLCHAR(XC,YC,1.2*CH,'Curvature vs X',0.0,16)
+ ENDIF
+C
+C--- Upper surface curvature
+ CALL GETCOLOR(ICOL0)
+ CALL NEWCOLORNAME('yellow')
+ CALL XYLINE(ILE,XB,CV,XOF,XSF,YOF,YSF,1)
+ XC = XSF*(XB(2*ILE/3)-XOF)
+ YC = YSF*(CV(2*ILE/3)-YOF)
+ CALL PLCHAR(XC+0.5*CH,YC+0.5*CH,CH,'Upper',0.0,5)
+C
+C--- LE curvature
+ CALL NEWCOLORNAME('red')
+ XC = XSF*(XLE-XOF)
+ YC = YSF*(CVLE-YOF)
+ CALL PLSYMB(XC,YC,CH,3,0.0,0)
+ CALL PLCHAR(XC+1.0*CH,YC-0.5*CH,CH,'LE',0.0,2)
+C
+C--- Lower surface curvature
+ CALL NEWCOLORNAME('cyan')
+ CALL XYLINE(NB-ILE+1,XB(ILE),CV(ILE),XOF,XSF,YOF,YSF,2)
+ XC = XSF*(XB((ILE+NB)/2)-XOF)
+ YC = YSF*(CV((ILE+NB)/2)-YOF)
+ CALL PLCHAR(XC+0.5*CH,YC+0.5*CH,CH,'Lower',0.0,5)
+C
+ CALL NEWCOLOR(ICOL0)
+ CALL PLFLUSH
+C
+C
+ 20 WRITE(*,*) ' '
+ WRITE(*,*) 'Airfoil curvature (yellow-upper, cyan-lower) '
+ IF(LCVEXP) THEN
+ WRITE(*,*) ' Range compressed using CV=(curvature)^n with n =',
+ & CVEXP
+ ENDIF
+ WRITE(*,*) ' CVLE = ',CVLE, ' at S = ',SBLE, ' at X = ',XLE
+ WRITE(*,*) ' CVmax = ',CVMAX,' at S = ',SVMAX,' at X = ',XVMAX
+ WRITE(*,*) ' CVmin = ',CVMIN,' at S = ',SVMIN,' at X = ',XVMIN
+C
+ WRITE(*,*) ' '
+ WRITE(*,*) 'Enter C for curvature plot'
+ WRITE(*,*) 'Enter N for curvature**N plot'
+ WRITE(*,*) 'Hit <cr> to exit'
+ ANSARG = ' '
+ CALL ASKC('..CPLO^',ANS,ANSARG)
+ IF(ANS.EQ.' ') RETURN
+C
+ RINPUT = 0.0
+ NINPUT = 1
+ CALL GETFLT(ANSARG,RINPUT,NINPUT,ERROR)
+C
+ IF(ANS.EQ.'n' .OR. ANS.EQ.'N') THEN
+ IF(NINPUT.GE.1) THEN
+ CVEXP = RINPUT
+ ELSE
+ CVEXP = 0.3
+ CALL ASKR('Enter curvature exponent (default 0.3)^',CVEXP)
+ ENDIF
+ LCVEXP = .TRUE.
+ GO TO 10
+ ENDIF
+ IF(ANS.EQ.'c' .OR. ANS.EQ.'C') THEN
+ LCVEXP = .FALSE.
+ GO TO 10
+ ENDIF
+ GO TO 20
+C
+ 1000 FORMAT(A)
+ END
+
+
+
+
+ SUBROUTINE CANG(X,Y,N,IPRINT, IMAX,AMAX)
+ DIMENSION X(*), Y(*)
+C-------------------------------------------------------------------
+C IPRINT=2: Displays all panel node corner angles
+C IPRINT=1: Displays max panel node corner angle
+C IPRINT=0: No display... just returns values
+C-------------------------------------------------------------------
+C
+ AMAX = 0.0
+ IMAX = 1
+C
+C---- go over each point, calculating corner angle
+ IF(IPRINT.EQ.2) WRITE(*,1050)
+ DO 30 I=2, N-1
+ DX1 = X(I) - X(I-1)
+ DY1 = Y(I) - Y(I-1)
+ DX2 = X(I) - X(I+1)
+ DY2 = Y(I) - Y(I+1)
+C
+C------ allow for doubled points
+ IF(DX1.EQ.0.0 .AND. DY1.EQ.0.0) THEN
+ DX1 = X(I) - X(I-2)
+ DY1 = Y(I) - Y(I-2)
+ ENDIF
+ IF(DX2.EQ.0.0 .AND. DY2.EQ.0.0) THEN
+ DX2 = X(I) - X(I+2)
+ DY2 = Y(I) - Y(I+2)
+ ENDIF
+C
+ CROSSP = (DX2*DY1 - DY2*DX1)
+ & / SQRT((DX1**2 + DY1**2) * (DX2**2 + DY2**2))
+ ANGL = ASIN(CROSSP)*(180.0/3.1415926)
+ IF(IPRINT.EQ.2) WRITE(*,1100) I, X(I), Y(I), ANGL
+ IF(ABS(ANGL) .GT. ABS(AMAX)) THEN
+ AMAX = ANGL
+ IMAX = I
+ ENDIF
+ 30 CONTINUE
+C
+ IF(IPRINT.GE.1) WRITE(*,1200) AMAX, IMAX, X(IMAX), Y(IMAX)
+C
+ RETURN
+C
+ 1050 FORMAT(/' i x y angle')
+CCC 120 0.2134 -0.0234 25.322
+ 1100 FORMAT(1X,I3, 2F9.4, F9.3)
+ 1200 FORMAT(/' Maximum panel corner angle =', F7.3,
+ & ' at i,x,y = ', I3, 2F9.4 )
+ END ! CANG
+
+
+
+ SUBROUTINE INTER(X0,XP0,Y0,YP0,S0,N0,SLE0,
+ & X1,XP1,Y1,YP1,S1,N1,SLE1,
+ & X,Y,N,FRAC)
+C .....................................................................
+C
+C Interpolates two source airfoil shapes into an "intermediate" shape.
+C
+C Procedure:
+C The interpolated x coordinate at a given normalized spline
+C parameter value is a weighted average of the two source
+C x coordinates at the same normalized spline parameter value.
+C Ditto for the y coordinates. The normalized spline parameter
+C runs from 0 at the leading edge to 1 at the trailing edge on
+C each surface.
+C .....................................................................
+C
+ REAL X0(N0),Y0(N0),XP0(N0),YP0(N0),S0(N0)
+ REAL X1(N1),Y1(N1),XP1(N1),YP1(N1),S1(N1)
+ REAL X(*),Y(*)
+C
+C---- number of points in interpolated airfoil is the same as in airfoil 0
+ N = N0
+C
+C---- interpolation weighting fractions
+ F0 = 1.0 - FRAC
+ F1 = FRAC
+C
+C---- top side spline parameter increments
+ TOPS0 = S0(1) - SLE0
+ TOPS1 = S1(1) - SLE1
+C
+C---- bottom side spline parameter increments
+ BOTS0 = S0(N0) - SLE0
+ BOTS1 = S1(N1) - SLE1
+C
+ DO 50 I=1, N
+C
+C------ normalized spline parameter is taken from airfoil 0 value
+ IF(S0(I).LT.SLE0) SN = (S0(I) - SLE0) / TOPS0 ! top side
+ IF(S0(I).GE.SLE0) SN = (S0(I) - SLE0) / BOTS0 ! bottom side
+C
+C------ set actual spline parameters
+ ST0 = S0(I)
+ IF(ST0.LT.SLE0) ST1 = SLE1 + TOPS1 * SN
+ IF(ST0.GE.SLE0) ST1 = SLE1 + BOTS1 * SN
+C
+C------ set input coordinates at common spline parameter location
+ XT0 = SEVAL(ST0,X0,XP0,S0,N0)
+ YT0 = SEVAL(ST0,Y0,YP0,S0,N0)
+ XT1 = SEVAL(ST1,X1,XP1,S1,N1)
+ YT1 = SEVAL(ST1,Y1,YP1,S1,N1)
+C
+C------ set interpolated x,y coordinates
+ X(I) = F0*XT0 + F1*XT1
+ Y(I) = F0*YT0 + F1*YT1
+C
+ 50 CONTINUE
+C
+ RETURN
+ END ! INTER
+
+
+
+ SUBROUTINE INTERX(X0,XP0,Y0,YP0,S0,N0,SLE0,
+ & X1,XP1,Y1,YP1,S1,N1,SLE1,
+ & X,Y,N,FRAC)
+C .....................................................................
+C
+C Interpolates two source airfoil shapes into an "intermediate" shape.
+C
+C Procedure:
+C The interpolated x coordinate at a given normalized spline
+C parameter value is a weighted average of the two source
+C x coordinates at the same normalized spline parameter value.
+C Ditto for the y coordinates. The normalized spline parameter
+C runs from 0 at the leading edge to 1 at the trailing edge on
+C each surface.
+C .....................................................................
+C
+ REAL X0(N0),Y0(N0),XP0(N0),YP0(N0),S0(N0)
+ REAL X1(N1),Y1(N1),XP1(N1),YP1(N1),S1(N1)
+ REAL X(N),Y(N)
+C
+C---- number of points in interpolated airfoil is the same as in airfoil 0
+ N = N0
+C
+C---- interpolation weighting fractions
+ F0 = 1.0 - FRAC
+ F1 = FRAC
+C
+ XLE0 = SEVAL(SLE0,X0,XP0,S0,N0)
+ XLE1 = SEVAL(SLE1,X1,XP1,S1,N1)
+C
+ DO 50 I=1, N
+C
+C------ normalized x parameter is taken from airfoil 0 value
+ IF(S0(I).LT.SLE0) XN = (X0(I) - XLE0) / (X0( 1) - XLE0)
+ IF(S0(I).GE.SLE0) XN = (X0(I) - XLE0) / (X0(N0) - XLE0)
+C
+C------ set target x and initial spline parameters
+ XT0 = X0(I)
+ ST0 = S0(I)
+ IF(ST0.LT.SLE0) THEN
+ XT1 = XLE1 + (X1( 1) - XLE1) * XN
+ ST1 = SLE1 + (S1( 1) - SLE1) * XN
+ ELSE
+ XT1 = XLE1 + (X1(N1) - XLE1) * XN
+ ST1 = SLE1 + (S1(N1) - SLE1) * XN
+ ENDIF
+C
+ CALL SINVRT(ST0,XT0,X0,XP0,S0,N0)
+ CALL SINVRT(ST1,XT1,X1,XP1,S1,N1)
+C
+C------ set input coordinates at common spline parameter location
+ XT0 = SEVAL(ST0,X0,XP0,S0,N0)
+ YT0 = SEVAL(ST0,Y0,YP0,S0,N0)
+ XT1 = SEVAL(ST1,X1,XP1,S1,N1)
+ YT1 = SEVAL(ST1,Y1,YP1,S1,N1)
+C
+C------ set interpolated x,y coordinates
+ X(I) = F0*XT0 + F1*XT1
+ Y(I) = F0*YT0 + F1*YT1
+C
+ 50 CONTINUE
+C
+ RETURN
+ END ! INTERX
+
+
+
+
+
+ SUBROUTINE BENDUMP(N,X,Y)
+ REAL X(*), Y(*)
+C
+ PEX = 16.0
+ CALL IJSECT(N,X,Y, PEX,
+ & AREA, SLEN,
+ & XMIN, XMAX, XEXINT,
+ & YMIN, YMAX, YEXINT,
+ & XC , YC ,
+ & XCT, YCT,
+ & AIXX , AIYY ,
+ & AIXXT, AIYYT,
+ & AJ , AJT )
+c CALL IJSECT(N,X,Y, PEX,
+c & AREA, SLEN,
+c & XC, XMIN, XMAX, XEXINT,
+c & YC, YMIN, YMAX, YEXINT,
+c & AIXX, AIXXT,
+c & AIYY, AIYYT,
+c & AJ , AJT )
+C
+ WRITE(*,*)
+ WRITE(*,1200) 'Area =', AREA
+ WRITE(*,1200) 'Slen =', SLEN
+ WRITE(*,*)
+ WRITE(*,1200) 'X-bending parameters(solid):'
+ WRITE(*,1200) ' Xc =', XC
+ WRITE(*,1200) ' max X-Xc =', XMAX-XC
+ WRITE(*,1200) ' min X-Xc =', XMIN-XC
+ WRITE(*,1200) ' Iyy =', AIYY
+ XBAR = MAX( ABS(XMAX-XC) , ABS(XMIN-XC) )
+ WRITE(*,1200) ' Iyy/(X-Xc)=', AIYY /XBAR
+ WRITE(*,*)
+ WRITE(*,1200) 'Y-bending parameters(solid):'
+ WRITE(*,1200) ' Yc =', YC
+ WRITE(*,1200) ' max Y-Yc =', YMAX-YC
+ WRITE(*,1200) ' min Y-Yc =', YMIN-YC
+ WRITE(*,1200) ' Ixx =', AIXX
+ YBAR = MAX( ABS(YMAX-YC) , ABS(YMIN-YC) )
+ WRITE(*,1200) ' Ixx/(Y-Yc)=', AIXX /YBAR
+ WRITE(*,*)
+ WRITE(*,1200) ' J =', AJ
+C
+ WRITE(*,*)
+ WRITE(*,*)
+ WRITE(*,1200) 'X-bending parameters(skin):'
+ WRITE(*,1200) ' Xc =', XCT
+ WRITE(*,1200) ' max X-Xc =', XMAX-XCT
+ WRITE(*,1200) ' min X-Xc =', XMIN-XCT
+ WRITE(*,1200) ' Iyy/t =', AIYYT
+ XBART = MAX( ABS(XMAX-XCT) , ABS(XMIN-XCT) )
+ WRITE(*,1200) ' Iyy/t(X-Xc)=', AIYYT /XBART
+ WRITE(*,*)
+ WRITE(*,1200) 'Y-bending parameters(skin):'
+ WRITE(*,1200) ' Yc =', YCT
+ WRITE(*,1200) ' max Y-Yc =', YMAX-YCT
+ WRITE(*,1200) ' min Y-Yc =', YMIN-YCT
+ WRITE(*,1200) ' Ixx/t =', AIXXT
+ YBART = MAX( ABS(YMAX-YCT) , ABS(YMIN-YCT) )
+ WRITE(*,1200) ' Ixx/t(Y-Yc)=', AIXXT /YBART
+ WRITE(*,*)
+ WRITE(*,1200) ' J/t =', AJT
+C
+c WRITE(*,*)
+c WRITE(*,1200) ' power-avg X-Xc =', XEXINT
+c WRITE(*,1200) ' power-avg Y-Yc =', YEXINT
+C
+ RETURN
+C
+ 1200 FORMAT(1X,A,G14.6)
+ END ! BENDUMP
+
+
+
+ SUBROUTINE BENDUMP2(N,X,Y,T)
+ REAL X(*), Y(*), T(*)
+C
+ DTR = ATAN(1.0) / 45.0
+C
+ PEX = 16.0
+ CALL IJSECT(N,X,Y, PEX,
+ & AREA, SLEN,
+ & XMIN, XMAX, XEXINT,
+ & YMIN, YMAX, YEXINT,
+ & XC , YC ,
+ & XCT, YCT,
+ & AIXX , AIYY ,
+ & AIXXT, AIYYT,
+ & AJ , AJT )
+c CALL IJSECT(N,X,Y, PEX,
+c & AREA, SLEN,
+c & XC, XMIN, XMAX, XEXINT,
+c & YC, YMIN, YMAX, YEXINT,
+c & AIXX, AIXXT,
+c & AIYY, AIYYT,
+c & AJ , AJT )
+C
+C
+ CALL AECALC(N,X,Y,T, 1,
+ & AREA,XCENA,YCENA,EI11A,EI22A,APX1A,APX2A)
+C
+ CALL AECALC(N,X,Y,T, 2,
+ & SLEN,XCENT,YCENT,EI11T,EI22T,APX1T,APX2T)
+C
+
+ WRITE(*,*)
+ WRITE(*,1200) 'Area =', AREA
+ WRITE(*,1200) 'Slen =', SLEN
+ WRITE(*,*)
+ WRITE(*,1200) 'X-bending parameters:'
+ WRITE(*,1200) 'solid centroid Xc=', XCENA
+ WRITE(*,1200) 'skin centroid Xc=', XCENT
+ WRITE(*,1200) ' solid max X-Xc =', XMAX-XCENA
+ WRITE(*,1200) ' solid min X-Xc =', XMIN-XCENA
+ WRITE(*,1200) ' skin max X-Xc =', XMAX-XCENT
+ WRITE(*,1200) ' skin min X-Xc =', XMIN-XCENT
+ WRITE(*,1200) ' solid Iyy =', EI22A
+ WRITE(*,1200) ' skin Iyy/t =', EI22T
+ XBARA = MAX( ABS(XMAX-XCENA) , ABS(XMIN-XCENA) )
+ XBART = MAX( ABS(XMAX-XCENT) , ABS(XMIN-XCENT) )
+ WRITE(*,1200) ' solid Iyy/(X-Xc)=', EI22A/XBARA
+ WRITE(*,1200) ' skin Iyy/t(X-Xc)=', EI22T/XBART
+C
+ WRITE(*,*)
+ WRITE(*,1200) 'Y-bending parameters:'
+ WRITE(*,1200) 'solid centroid Yc=', YCENA
+ WRITE(*,1200) 'skin centroid Yc=', YCENT
+ WRITE(*,1200) ' solid max Y-Yc =', YMAX-YCENA
+ WRITE(*,1200) ' solid min Y-Yc =', YMIN-YCENA
+ WRITE(*,1200) ' skin max Y-Yc =', YMAX-YCENT
+ WRITE(*,1200) ' skin min Y-Yc =', YMIN-YCENT
+ WRITE(*,1200) ' solid Ixx =', EI11A
+ WRITE(*,1200) ' skin Ixx/t =', EI11T
+ YBARA = MAX( ABS(YMAX-YCENA) , ABS(YMIN-YCENA) )
+ YBART = MAX( ABS(YMAX-YCENT) , ABS(YMIN-YCENT) )
+ WRITE(*,1200) ' solid Ixx/(Y-Yc)=', EI11A/YBARA
+ WRITE(*,1200) ' skin Ixx/t(Y-Yc)=', EI11T/YBART
+C
+ WRITE(*,*)
+ WRITE(*,1200) ' solid principal axis angle (deg ccw) =', APX1A/DTR
+ WRITE(*,1200) ' skin principal axis angle (deg ccw) =', APX1T/DTR
+
+c WRITE(*,*)
+c WRITE(*,1200) ' power-avg X-Xc =', XEXINT
+c WRITE(*,1200) ' power-avg Y-Yc =', YEXINT
+C
+ WRITE(*,*)
+ WRITE(*,1200) ' solid J =', AJ
+ WRITE(*,1200) ' skin J/t =', AJT
+ RETURN
+C
+ 1200 FORMAT(1X,A,G14.6)
+ END ! BENDUMP2
+
+
+
+ SUBROUTINE IJSECT(N,X,Y, PEX,
+ & AREA, SLEN,
+ & XMIN, XMAX, XEXINT,
+ & YMIN, YMAX, YEXINT,
+ & XC , YC ,
+ & XCT, YCT,
+ & AIXX , AIYY ,
+ & AIXXT, AIYYT,
+ & AJ , AJT )
+ DIMENSION X(*), Y(*)
+C
+ XMIN = X(1)
+ XMAX = X(1)
+ YMIN = Y(1)
+ YMAX = Y(1)
+C
+ DX = X(1) - X(N)
+ DY = Y(1) - Y(N)
+ DS = SQRT(DX*DX + DY*DY)
+ XAVG = 0.5*(X(1) + X(N))
+ YAVG = 0.5*(Y(1) + Y(N))
+C
+ X_DY = DY * XAVG
+ XX_DY = DY * XAVG**2
+ XXX_DY = DY * XAVG**3
+ X_DS = DS * XAVG
+ XX_DS = DS * XAVG**2
+C
+ Y_DX = DX * YAVG
+ YY_DX = DX * YAVG**2
+ YYY_DX = DX * YAVG**3
+ Y_DS = DS * YAVG
+ YY_DS = DS * YAVG**2
+C
+ C_DS = DS
+C
+ DO 10 I = 2, N
+ DX = X(I) - X(I-1)
+ DY = Y(I) - Y(I-1)
+ DS = SQRT(DX*DX + DY*DY)
+ XAVG = 0.5*(X(I) + X(I-1))
+ YAVG = 0.5*(Y(I) + Y(I-1))
+C
+ X_DY = X_DY + DY * XAVG
+ XX_DY = XX_DY + DY * XAVG**2
+ XXX_DY = XXX_DY + DY * XAVG**3
+ X_DS = X_DS + DS * XAVG
+ XX_DS = XX_DS + DS * XAVG**2
+C
+ Y_DX = Y_DX + DX * YAVG
+ YY_DX = YY_DX + DX * YAVG**2
+ YYY_DX = YYY_DX + DX * YAVG**3
+ Y_DS = Y_DS + DS * YAVG
+ YY_DS = YY_DS + DS * YAVG**2
+C
+ C_DS = C_DS + DS
+C
+ XMIN = MIN(XMIN,X(I))
+ XMAX = MAX(XMAX,X(I))
+ YMIN = MIN(YMIN,Y(I))
+ YMAX = MAX(YMAX,Y(I))
+ 10 CONTINUE
+C
+ AREA = -Y_DX
+ SLEN = C_DS
+C
+ IF(AREA.EQ.0.0) RETURN
+C
+ XC = XX_DY / (2.0*X_DY)
+ XCT = X_DS / C_DS
+ AIYY = XXX_DY/3.0 - XX_DY*XC + X_DY*XC**2
+ AIYYT = XX_DS - X_DS*XCT*2.0 + C_DS*XCT**2
+C
+ YC = YY_DX / (2.0*Y_DX)
+ YCT = Y_DS / C_DS
+ AIXX = -YYY_DX/3.0 + YY_DX*YC - Y_DX*YC**2
+ AIXXT = YY_DS - Y_DS*YCT*2.0 + C_DS*YCT**2
+C
+C
+ SINT = 0.
+ XINT = 0.
+ YINT = 0.
+C
+ DO 20 I=2, N
+ DX = X(I) - X(I-1)
+ DY = Y(I) - Y(I-1)
+ DS = SQRT(DX*DX + DY*DY)
+ XAVG = 0.5*(X(I) + X(I-1)) - XC
+ YAVG = 0.5*(Y(I) + Y(I-1)) - YC
+C
+ SINT = SINT + DS
+cc XINT = XINT + DS * ABS(XAVG)**PEX
+cc YINT = YINT + DS * ABS(YAVG)**PEX
+ 20 CONTINUE
+C
+ DO I=1, N-1
+ IF(X(I+1) .GE. X(I)) GO TO 30
+ ENDDO
+ IMID = N/2
+ 30 IMID = I
+C
+ AJ = 0.0
+ DO I = 2, IMID
+ XAVG = 0.5*(X(I) + X(I-1))
+ YAVG = 0.5*(Y(I) + Y(I-1))
+ DX = X(I-1) - X(I)
+C
+ IF(XAVG.GT.X(N)) THEN
+ YOPP = Y(N)
+ GO TO 41
+ ENDIF
+ IF(XAVG.LE.X(IMID)) THEN
+ YOPP = Y(IMID)
+ GO TO 41
+ ENDIF
+C
+ DO J = N, IMID, -1
+ IF(XAVG.GT.X(J-1) .AND. XAVG.LE.X(J)) THEN
+ FRAC = (XAVG - X(J-1))
+ & / (X(J) - X(J-1))
+ YOPP = Y(J-1) + (Y(J)-Y(J-1))*FRAC
+ GO TO 41
+ ENDIF
+ ENDDO
+ 41 CONTINUE
+C
+ AJ = AJ + ABS(YAVG-YOPP)**3 * DX / 3.0
+ ENDDO
+C
+ AJT = 4.0*AREA**2/SLEN
+C
+cc XEXINT = (XINT/SINT)**(1.0/PEX)
+cc YEXINT = (YINT/SINT)**(1.0/PEX)
+C
+ RETURN
+ END ! IJSECT
+C
+
+
+ SUBROUTINE AREFINE(X,Y,S,XS,YS,N, ATOL,
+ & NDIM,NNEW,XNEW,YNEW,X1,X2)
+C-------------------------------------------------------------
+C Adds points to a x,y spline contour wherever
+C the angle between adjacent segments at a node
+C exceeds a specified threshold. The points are
+C added 1/3 of a segment before and after the
+C offending node.
+C
+C The point adding is done only within X1..X2.
+C
+C Intended for doubling the number of points
+C of Eppler and Selig airfoils so that they are
+C suitable for clean interpolation using Xfoil's
+C arc-length spline routines.
+C------------------------------------------------------
+ REAL X(*), Y(*), S(*), XS(*), YS(*)
+ REAL XNEW(NDIM), YNEW(NDIM)
+ LOGICAL LREF
+C
+ ATOLR = ATOL * 3.14159/180.0
+C
+ K = 1
+ XNEW(K) = X(1)
+ YNEW(K) = Y(1)
+C
+ DO 10 I = 2, N-1
+ IM = I-1
+ IP = I+1
+C
+ DXM = X(I) - X(I-1)
+ DYM = Y(I) - Y(I-1)
+ DXP = X(I+1) - X(I)
+ DYP = Y(I+1) - Y(I)
+C
+ CRSP = DXM*DYP - DYM*DXP
+ DOTP = DXM*DXP + DYM*DYP
+ IF(CRSP.EQ.0.0 .AND. DOTP.EQ.0.0) THEN
+ ASEG = 0.0
+ ELSE
+ ASEG = ATAN2( CRSP , DOTP )
+ ENDIF
+C
+ LREF = ABS(ASEG) .GT. ATOLR
+C
+ IF(LREF) THEN
+C------- add extra point just before this node
+ SMID = S(I) - 0.3333*(S(I)-S(I-1))
+ XK = SEVAL(SMID,X,XS,S,N)
+ YK = SEVAL(SMID,Y,YS,S,N)
+ IF(XK.GE.X1 .AND. XK.LE.X2) THEN
+ K = K + 1
+ IF(K .GT. NDIM) GO TO 90
+ XNEW(K) = XK
+ YNEW(K) = YK
+ ENDIF
+ ENDIF
+C
+C------ add the node itself
+ K = K + 1
+ IF(K .GT. NDIM) GO TO 90
+ XNEW(K) = X(I)
+ YNEW(K) = Y(I)
+C
+ IF(LREF) THEN
+C------- add extra point just after this node
+ SMID = S(I) + 0.3333*(S(I+1)-S(I))
+ XK = SEVAL(SMID,X,XS,S,N)
+ YK = SEVAL(SMID,Y,YS,S,N)
+ IF(XK.GE.X1 .AND. XK.LE.X2) THEN
+ K = K + 1
+ IF(K .GT. NDIM) GO TO 90
+ XNEW(K) = XK
+ YNEW(K) = YK
+ ENDIF
+ ENDIF
+ 10 CONTINUE
+C
+ K = K + 1
+ IF(K .GT. NDIM) GO TO 90
+ XNEW(K) = X(N)
+ YNEW(K) = Y(N)
+C
+ NNEW = K
+ RETURN
+C
+ 90 CONTINUE
+ WRITE(*,*) 'SDOUBLE: Arrays will overflow. No action taken.'
+ NNEW = 0
+ RETURN
+C
+ END ! AREFINE
+
+
+ SUBROUTINE SCHECK(X,Y,N, STOL, LCHANGE)
+C-------------------------------------------------------------
+C Removes points from an x,y spline contour wherever
+C the size of a segment between nodes falls below a
+C a specified threshold of the adjacent segments.
+C The two node points defining the short segment are
+C replaced with a single node at their midpoint.
+C Note that the number of nodes may be altered by
+C this routine.
+C
+C Intended for eliminating odd "micro" panels
+C that occur when blending a flap to a foil.
+C If LCHANGE is set on return the airfoil definition
+C has been changed and resplining should be done.
+C
+C The recommended value for STOL is 0.05 (meaning
+C segments less than 5% of the length of either adjoining
+C segment are removed). 4/24/01 HHY
+C------------------------------------------------------
+ REAL X(*), Y(*)
+ LOGICAL LCHANGE
+C
+ LCHANGE = .FALSE.
+C--- Check STOL for sanity
+ IF(STOL.GT.0.3) THEN
+ WRITE(*,*) 'SCHECK: Bad value for small panels (STOL > 0.3)'
+ RETURN
+ ENDIF
+C
+ 10 DO 20 I = 2, N-2
+ IM1 = I-1
+ IP1 = I+1
+ IP2 = I+2
+C
+ DXM1 = X(I) - X(I-1)
+ DYM1 = Y(I) - Y(I-1)
+ DSM1 = SQRT(DXM1*DXM1 + DYM1*DYM1)
+C
+ DXP1 = X(I+1) - X(I)
+ DYP1 = Y(I+1) - Y(I)
+ DSP1 = SQRT(DXP1*DXP1 + DYP1*DYP1)
+C
+ DXP2 = X(I+2) - X(I+1)
+ DYP2 = Y(I+2) - Y(I+1)
+ DSP2 = SQRT(DXP2*DXP2 + DYP2*DYP2)
+C
+C------- Don't mess with doubled points (slope breaks)
+ IF(DSP1.EQ.0.0) GO TO 20
+C
+ IF(DSP1.LT.STOL*DSM1 .OR. DSP1.LT.STOL*DSP2) THEN
+C------- Replace node I with average of I and I+1
+ X(I) = 0.5*(X(I)+X(I+1))
+ Y(I) = 0.5*(Y(I)+Y(I+1))
+C------- Remove node I+1
+ DO L = I+1, N
+ X(L) = X(L+1)
+ Y(L) = Y(L+1)
+ END DO
+ N = N - 1
+ LCHANGE = .TRUE.
+ WRITE(*,*) 'SCHECK segment removed at ',I
+ GO TO 10
+ ENDIF
+C
+ 20 CONTINUE
+C
+ RETURN
+ END ! SCHECK
+
+
+
+ SUBROUTINE HALF(X,Y,S,N)
+C-------------------------------------------------
+C Halves the number of points in airfoil
+C-------------------------------------------------
+ REAL X(*), Y(*), S(*)
+C
+ K = 1
+ INEXT = 3
+ DO 20 I=2, N-1
+C------ if corner is found, preserve it.
+ IF(S(I) .EQ. S(I+1)) THEN
+ K = K+1
+ X(K) = X(I)
+ Y(K) = Y(I)
+ K = K+1
+ X(K) = X(I+1)
+ Y(K) = Y(I+1)
+ INEXT = I+3
+ ENDIF
+C
+ IF(I.EQ.INEXT) THEN
+ K = K+1
+ X(K) = X(I)
+ Y(K) = Y(I)
+ INEXT = I+2
+ ENDIF
+C
+ 20 CONTINUE
+ K = K+1
+ X(K) = X(N)
+ Y(K) = Y(N)
+C
+C---- set new number of points
+ N = K
+C
+ RETURN
+ END ! HALF
+
+
+