subroutine LVFoil(XB,YB,n, alphaD, X, Y, Cp)
C
C Linear vorticity surface panel method for airfoils.
C
C Adapted from Kuethe and Chow 4th Edition
C This version by I. Kroo 2-2-87
C
C No attempt has been made to make this particularly fast or efficient.
C
C
C Inputs:
C -------
C XB, YB x and y coordinates of panel edges starting at trailing edge
C proceeding forward on lower surface, wrapping around
C leading edge and then running back to the trailingedge.
C n Number of panels
C alphaD Angle of attack in degrees
C
C Outputs:
C --------
C X, Y coordinates of panel centers
C Cp incompressible pressure coefficient at panel center
C
C Declarations:
C -------------
implicit none
integer ndim
parameter (ndim = 100)
real XB(ndim), YB(ndim), X(ndim), Y(ndim), S(ndim), SINT(ndim)
real COST(ndim), THETA(ndim), V(ndim), RHS(ndim), CP(ndim)
real CN1(ndim,ndim), CN2(ndim,ndim)
real CT1(ndim,ndim), CT2(ndim,ndim)
real AN(ndim,ndim), AT(ndim,ndim)
real pi, alpha, A, B, C, D, E, F, G, P, Q, P1, P2, P3, P4
real COSA, SINA, ANmax, alphaD
integer n, np1, i, j, IP(ndim), ip1
logical err
C
C Constants
C ---------
pi = 3.14159265
C
C Calculation of geometric data:
C ------------------------------
np1 = n+1
do 1 i = 1,n
ip1 = i+1
X(i) = .5 * (XB(i) + XB(ip1))
Y(i) = .5 * (YB(i) + YB(ip1))
S(i) = SQRT( (XB(ip1)-XB(i))**2 + (YB(ip1)-YB(i))**2)
THETA(i) = ATAN2( (YB(ip1)-YB(i)), (XB(ip1)-XB(i)) )
SINT(i) = SIN(THETA(i))
COST(i) = COS(THETA(i))
1 continue
alpha = alphaD *pi/180.
SINA = SIN(alpha)
COSA = COS(alpha)
C
C Calculation of influence coefficients
C -------------------------------------
C
do 3 i = 1,n
do 2 j = 1,n
if (i.eq.j) then
CN1(i,j) = -1.
CN2(i,j) = 1.
CT1(i,j) = .5*pi
CT2(i,j) = CT1(i,j)
else
A = -(X(i)-XB(j))*COST(j) - (Y(i)-YB(j))*SINT(j)
B = (X(i)-XB(j))**2 + (Y(i)-YB(j))**2
C = SINT(i)*COST(j)-COST(i)*SINT(j)
D = COST(i)*COST(j)+SINT(i)*SINT(j)
E = (X(i)-XB(j))*SINT(j) - (Y(i)-YB(j))*COST(j)
F = ALOG( 1. + S(j)*(S(j)+2.*A)/B )
G = ATAN2(E*S(j), B+A*S(j))
P1 = 1.-2.*SINT(j)*SINT(j)
P2 = 2.*SINT(j)*COST(j)
P3 = SINT(i)*P1 - COST(i)*P2
P4 = COST(i)*P1 + SINT(i)*P2
P = (X(i)-XB(j)) * P3 + (Y(i)-YB(j)) * P4
Q = (X(i)-XB(j)) * P4 - (Y(i)-YB(j)) * P3
CN2(i,j) = D + ( .5*Q*F - (A*C+D*E)*G )/S(j)
CN1(i,j) = .5*D*F + C*G - CN2(i,j)
CT2(i,j) = C + ( .5*P*F + (A*D-C*E)*G )/S(j)
CT1(i,j) = .5*C*F - D*G - CT2(i,j)
end if
2 continue
3 continue
ANmax = 0.
do 5 i = 1,n
AN(i,1) = CN1(i,1)
AN(i,np1) = CN2(i,n)
AT(i,1) = CT1(i,1)
AT(i,np1) = CT2(i,n)
do 4 j = 2,n
AN(i,j) = CN1(i,j) + CN2(i,j-1)
AT(i,j) = CT1(i,j) + CT2(i,j-1)
if(ABS(AN(i,j)).gt.ANmax)ANmax = ABS(AN(i,j))
4 continue
5 continue
C The Kutta Condition is imposed by the relation:
C A*gamma(1) + A*gamma(n+1) = 0. A would be 1 but for matrix
C conditioning problems.
AN(np1,1) = 1
AN(np1,np1) = 1
do 6 j = 2,n
AN(np1,j) = 0.
6 continue
RHS(np1) = 0.
C
C Decompose and solve the system:
C -------------------------------
C
call DECOMP(np1, ndim, AN, IP)
C RHS(i) = SIN( THETA(i)-alpha )
do 11 i = 1,n
RHS(i) = SINT(i)*COSA - COST(i)*SINA
11 continue
RHS(np1) = 0.
call SOLVE(np1, ndim, AN, RHS, IP)
C
C Compute derived quantities:
C ---------------------------
C
C V(i) = COS( THETA(i) - alpha )
do 8 i = 1,n
V(i) = COST(i)*COSA + SINT(i)*SINA
do 7 j = 1,np1
V(i) = V(i) + AT(i,j) * RHS(j)
7 continue
Cp(i) = 1.-V(i)*V(i)
8 continue
return
end
C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION
C
C
C Variables:
C N The size of the system to be solved
C NDIM The dimension of the maxtrix as declared in the calling program
C C(NDIM,NDIM) The matrix to be decomposed by L-U decomposition.
C On return C is the decomposed form of the matrix.
C IP(NDIM) A vector of pivots which is used by the routine SOLVE.
C If IP(1)=0 then the matrix is singular.
C
SUBROUTINE DECOMP(N,NDIM,C,IP)
REAL C(NDIM,NDIM),T
INTEGER IP(NDIM)
IP(N)=1
DO 6 K=1,N
IF(K.EQ.N)GOTO 5
KP1=K+1
M=K
DO 1 I=KP1,N
IF(ABS(C(I,K)).GT. ABS(C(M,K)))M=I
1 CONTINUE
IP(K)=M
IF(M.NE.K)IP(N)=-IP(N)
T=C(M,K)
C(M,K)=C(K,K)
C(K,K)=T
IF (T.EQ.0.0)GOTO 5
DO 2 I=KP1,N
C(I,K)=-C(I,K)/T
2 CONTINUE
DO 4 J=KP1,N
T=C(M,J)
C(M,J)=C(K,J)
C(K,J)=T
IF (T.EQ.0.0)GOTO 4
DO 3 I=KP1,N
C(I,J)=C(I,J)+C(I,K)*T
3 CONTINUE
4 CONTINUE
5 IF(C(K,K).EQ.0.0)IP(N)=0
6 CONTINUE
RETURN
END
C SOLUTION OF LINEAR SYSTEM CX = B
C
C Variables:
C N The size of the system to be solved
C NDIM The dimension of the maxtrix as declared in the calling program
C C(NDIM,NDIM) The matrix as returned by the routine DECOMP
C IP(NDIM) A vector of pivots. If IP(1)=0 then the matrix is singular.
C B(NDIM) The right hand side vector. On return, the solution vector.
C
C
SUBROUTINE SOLVE(N,NDIM,C,B,IP)
REAL C(NDIM,NDIM),B(NDIM),T
INTEGER IP(NDIM)
IF(N.NE.1)THEN
NM1=N-1
DO 2 K = 1,NM1
KP1=K+1
M=IP(K)
T=B(M)
B(M)=B(K)
B(K)=T
DO 1 I=KP1,N
B(I)=B(I)+C(I,K)*T
1 CONTINUE
2 CONTINUE
DO 4 KB=1,NM1
KM1=N-KB
K=KM1+1
B(K)=B(K)/C(K,K)
T=-B(K)
DO 3 I=1,KM1
B(I)=B(I)+C(I,K)*T
3 CONTINUE
4 CONTINUE
END IF
B(1)=B(1)/C(1,1)
RETURN
END