C      REVISED 12/12/85   JPS
C  AERODYNAMIC PROPELLER ANALYSIS
C  METHOD- STRIP THEORY WITH GOLDSTEIN INDUCTION FACTOR
C  SEE NOTES FOR DERIVATION AND EQUATIONS
C
C  INPUT FORMAT
C  CARD #1
C  M,IPRINT       (FORMAT(3I10)
C  CARD#2,M+1
C  X,COR,BETA,TOC,CD   (FORMAT(7F13.6))
C  CARD #M+2
C  RPM,B,VIN,DIA,RHO     (FORMAT(7F13.6))
C
C
C  VARIABLE DEFINITIONS
C       M - NUMBER OF RADIAL STATIONS
C       IPRINT - PRINT PARAMETER ( =0,FULL OUTPUT ;=1, SUMMARY OUTPUT ONLY)
C       X - DIMENSIONLESS RADIAL POSITION,R/R,(MUST BE .3,.35,.4,.45--.95,1.)
C       COR - CHORD OVER RADIUS, C/R
C       BETA - TWIST DISTRIBUTION- ANGLE PLANE OF ROTATION TO ZERO LIFT LINE
C       TOC - THICKNESS OVER CHORD, T/C
C       CD - DRAG COEFFCIENT
C       RPM - REVOLUTIONS PER MINUTE
C       B - NUMBER OF BLADES
C       VIN - VELOCITY AT INFINITY ( FLIGHT SPEED ) (FEET/SECOND)
C       DIA - DIAMETER (inches)
C       RHO - DENSITY OF AIR (SLUGS PER CUBIC FOOT) (.002378 AT STP)
C	BETAIN- BETA 3/4 OF INPUT TWIST DISTRIBUTION
C	BETOUT- BETA 3/4 DESIRED FOR CALCULATION
C
C
        PROGRAM GOLD
        DIMENSION DCTDX(20),DCPDX(20),GUESS(20),ERROR(20)
        DIMENSION X(20),COR(20),BETA(20),TOC(20),CD(20),CL(20)
	DIMENSION BETAI(20),HEADER(20)
	character*64 filnam,fnam2
	character IY
        REAL J,KAPPA,KAPPA2,KAPPA3,KAPPA4
        INTEGER OUTPUT
        DATA INPUT,OUTPUT,PI,DEG2RA/12,11,3.14159,.017453/
C
	open(11,file='gold.out')
c
C   GET NAME AND OPEN FILE FOR INPUT
C
        WRITE(*,5000)
5000    FORMAT(2X,'FILE NAME OF INPUT DATA ?  ',/)
	READ(*,5001)FILNAM
5001	FORMAT(A)
        open(INPUT,file=FILNAM)
	WRITE(OUTPUT,5002)FILNAM
5002	FORMAT(5X,'NAME OF INPUT DATA FILE -  ',A)
	WRITE(*,5800)
5800	FORMAT(5X,' INPUT COMMENTS FOR HEADER OF PRINTED OUTPUT
     1  -80 CHARS MAX')
	READ(*,5801)(HEADER(I), I=1,20)
5801	FORMAT(20A4)
	WRITE(OUTPUT,5802)(HEADER(I), I=1,20)
5802	FORMAT(5X,20A4)
C
C   GET NAME AND OPEN DISC FILE FOR OUTPUT
C
	WRITE(*,5010)
5010	FORMAT(5X,'DO YOU WANT OUTPUT WRITTEN TO DISC ? (y OR n)')
	READ(*,5011)IY
5011	FORMAT(A)
	IF(IY.Ne.'y')gO TO 5012
	WRITE(*,5013)
5013	FORMAT(5X,'INPUT FILE NAME FOR OUTPUT',/)
	IFILE=13
	read(*,5001)fnam2
	open(IFILE,file=fnam2,status='new')
5012	CONTINUE
C
C
        READ (INPUT,1000)M,IPRINT
1000    FORMAT(2I10)
        WRITE(OUTPUT,6000)
6000    FORMAT(20X,'INPUT DATA',//)
        WRITE(OUTPUT,2000)M,IPRINT
2000    FORMAT(5X,'NUMBER OF RADIAL STATIONS = ',I5,'  PRINT PARAMETER
     1  =  ',I5)
	WRITE(OUTPUT,5003)
5003	FORMAT(8X,'X',10X,'COR',11X,'BETA',12X,'TOC',10X,'CD')
        DO 10 I=1,M
        READ(INPUT,1100)X(I),COR(I),BETAI(I),TOC(I),CD(I)
        WRITE(OUTPUT,2100)X(I),COR(I),BETAI(I),TOC(I),CD(I)
10      CONTINUE
1100    FORMAT(7F13.6)
2100    FORMAT(1X,7F13.6)
        IF(IPRINT.EQ.1)WRITE(OUTPUT,6200)
20      READ(INPUT,1100,END=500)RPM,B,VIN,DIA,RHO,BETAIN,BETOUT
C
C   ROTATE BLADE TO DESIRED BETA 3/4 (BETOUT) FROM
C   THE BETA 3/4 OF THE INPUT TWIST DISTRIBUTION (BETAIN)
C
	DO 23 I=1,M
23	BETA(I)=( BETAI(I) +BETOUT - BETAIN)*DEG2RA
C
C
	IB=IFIX(B+.5)
        RPS=RPM/60.
        J=VIN/RPS/DIA*12.
        VINWR=J/PI
        IF (IPRINT.EQ.1)GO TO 21
	WRITE(OUTPUT,5803)(HEADER(I), I=1,20)
5803	FORMAT('1',20A4)
	WRITE(OUTPUT,5001)FILNAM
	WRITE(OUTPUT,5004)
5004	FORMAT(8X,'J',7X,'# BLADE',7X,'VELOCITY',7X,'DIA',9X,
     1  'DENSITY',7X,'BETAIN',10X,'BETOUT')
        WRITE(OUTPUT,2100)J,B,VIN,DIA,RHO,BETAIN,BETOUT
        WRITE(OUTPUT,6100)
6100    FORMAT(//,5X,'X      BETA    PHI     KAPPA    ALPHAI',9X,
     1  'WT/WR',11X,'WA/WR',8X,'VE/WR',10X,'CL',10X,'DCTDX',10X,'DCPDX'
     2  ,/)
21      CONTINUE
C
C   BEGIN CALCULATION
C
        DO 30 I=1,M
        PHI=ATAN(VINWR/X(I))
C
C  SEE FIG. 4 FOR ANGLE DEFINITIONS
C  SIGMA - SOLIDITY
C
        SIGMA=B*COR(I)/PI
C
C   LAST STATION AT TIP
	IF (I.EQ.15) GO TO 65
C  BEGIN ITERATION
        ITER=1
C
C  FIRST GUESS SET TO AVOID SQRT OF NEGATIVE NUMBER
	GUESS(ITER)=(X(I)-SQRT(X(I)**2+VINWR**2))/2 + .00001
40      WTWR=GUESS(ITER)
C  EQUATION #4
c	 print *,X(I),WTWR,ITER,kappa
        ARG=VINWR**2+4.*WTWR*(X(I)-WTWR)
	IF(ARG)95,95,96
95	WAWR=WTWR*10
	GO TO 97
96      WAWR=.5*(-VINWR+SQRT(VINWR**2+4.*WTWR*(X(I)-WTWR)))
97	CONTINUE
C  EQUATION #5
        VEWR=SQRT((VINWR+WAWR)**2+(X(I)-WTWR)**2)
C
        ALPHAI=ATAN(WTWR/WAWR)-PHI
C  EQUATION #3
        CL(I)=2*PI*(BETA(I)-ALPHAI-PHI)
        PHIALP=PHI+ALPHAI
C OBTAIN GOLDSTEIN KAPPA FACTOR FROM SUBFUNCTION KAPPA2,KAPPA3,KAPPA4
	GO TO (41,42,43) IB-1
41      KAPPA=KAPPA2(X(I),PHIALP)
c	print *,kappa
	GO TO 44
42	KAPPA=KAPPA3(X(I),PHIALP)
	GO TO 44
43	KAPPA=KAPPA4(X(I),PHIALP)
44	CONTINUE
C  EQUATION #2
        ERROR(ITER)=SIGMA*CL(I)*VEWR-8.*X(I)*KAPPA*WTWR
        IF(ABS(ERROR(ITER))-.00001)60,60,50
50      ITER=ITER+1
        IF(ITER.GE.20)GO TO 60
C NO MORE THAN 20 ITERATIONS ( NORMALLY 4-6)
        IF(ITER.GT.2)GO TO 55
C SECOND GUESS ONLY
        GUESS(ITER)=GUESS(ITER-1)+.01
        GO TO 40
55      CONTINUE
C  NEWTON'S METHOD FOR GUESS> 2
	SLOPE=(ERROR(ITER-2)-ERROR(ITER-1))/(GUESS(ITER-2)-GUESS(ITER-1))
C
        BINTER=ERROR(ITER-2)-SLOPE*GUESS(ITER-2)
        GUESS(ITER)=-BINTER/SLOPE
        GO TO 40
C
C   AT X=1.0 CL=0.0
C
65	CL(I)=0.0
	ALPHAI=BETA(I)-PHI
	VRWR=SQRT(VINWR**2+1)
	WAWR=VRWR*SIN(ALPHAI)*COS(ALPHAI+PHI)
	WTWR=VRWR*SIN(ALPHAI)*SIN(PHI+ALPHAI)
        VEWR=SQRT((VINWR+WAWR)**2+(X(I)-WTWR)**2)
	KAPPA=0.0
60      CONTINUE
        WRITE(*,2200)ITER
2200    FORMAT(1X,'THE NUMBER OF ITERATIONS =  ',I10)
        C1=PI*SIGMA*(J**2+(PI*X(I))**2)/8.
        COS1=COS(PHI+ALPHAI)
        SIN1=SIN(PHI+ALPHAI)
C  EQUATION #1A  (INTEGRAND ONLY)
        DCTDX(I)=C1*(CL(I)*COS1-CD(I)*SIN1)
C  EQUATION #1B  (INTEGRAND ONLY)
        DCPDX(I)=C1*PI*X(I)*(CL(I)*SIN1+CD(I)*COS1)
        IF (IPRINT.EQ.1)GO TO 30
C  OUTPUT RADIAL DISTRIBUTION INFORMATION (SKIP IF IPRINT=1)
        WRITE(OUTPUT,2300)X(I),BETA(I),PHI,KAPPA,ALPHAI,WTWR,WAWR,VEWR
     1  ,CL(I),DCTDX(I),DCPDX(I)
C  WRITE OUTPUT TO DISC  IF FILE OPENED
	IF(IY.NE.'y')GO TO 30
        WRITE(IFILE,2350)X(I),BETA(I),PHI,KAPPA,ALPHAI,WTWR,WAWR,VEWR
     1  ,CL(I),DCTDX(I),DCPDX(I)
2300    FORMAT(1X,4(2X,F6.3),7(3X,E11.4))
2350	FORMAT(12E11.4)
30      CONTINUE
C
C  END OF MAIN CALUCLATION LOOP
C
        IF(IPRINT.EQ.0)WRITE(OUTPUT,6200)
6200    FORMAT(//,7X,'J',9X,'CT',9X,'CP',7X,'EFF',4X,'VIN',8X,'RPM',
     1  7X,'THRUST',4X,'POWER',4X,'# BLADES',4X,'BETA 3/4',2X,'THRUST/
     2  POWER',/)
C
C  INTERGRATE USING TRAPEZOID RULE FOR TOTL POWER AND THRUS
C
        CTSUM=0.0
        CPSUM=0.0
        DO 70 I=2,M-1
        CPSUM=CPSUM+DCPDX(I)
        CTSUM=CTSUM+DCTDX(I)
70      CONTINUE
        STEP=X(2)-X(1)
        CP=STEP*(CPSUM+(DCPDX(1)+DCPDX(M))/2.)
        CT=STEP*(CTSUM+(DCTDX(1)+DCTDX(M))/2.)
        ETA=CT*J/CP
        C10=RHO*(RPS)**2*(DIA/12.)**4
        THRUST=CT*C10
        POWER=CP*C10*RPS*DIA/550./12.
	IF (POWER)73,72,73
72	TOHP=0.0
	GO TO 74
73	TOHP=THRUST/POWER
74	CONTINUE
	write(20,2401)j
	write(21,2401)ct
	write(22,2401)cp
	write(23,2401)eta
2401    format(f10.6)
        WRITE(OUTPUT,2400)J,CT,CP,ETA,VIN,RPM,THRUST,POWER,B,BETOUT,TOHP
	IF (IY.NE.'y')GO TO 71
        WRITE(IFILE,2350)J,CT,CP,ETA,VIN,RPM,THRUST,POWER,B,BETOUT,TOHP
2400    FORMAT(1X,12F10.4,/)
71	CONTINUE
        GO TO 20
500     STOP
        END
        FUNCTION KAPPA2(X,PHIALP)
C
C   RETURNS THE GOLDSTEIN KAPPA FACTOR FOR TWO BLADED PROPELLERS
C
        REAL KAPPA2
        DIMENSION F(11,14)
        DATA F/
     1  1.00,1.00,.994,.978,.958,.944,.930,.922,.916,.935,1.02,
     2  1.00,1.00,.993,.961,.940,.912,.880,.855,.837,.843,.891,
     3  1.00,.999,.992,.955,.923,.880,.833,.788,.757,.750,.762,
     4  1.00,.998,.991,.959,.906,.848,.784,.722,.677,.657,.633,
     5  1.00,.998,.981,.930,.865,.794,.725,.663,.618,.592,.563,
     6  1.00,.997,.971,.902,.824,.742,.666,.605,.560,.528,.494,
     7  1.00,.997,.961,.874,.783,.690,.608,.547,.502,.464,.425,
     8  .999,.992,.931,.824,.723,.627,.550,.494,.448,.412,.375,
     9  .999,.988,.901,.774,.663,.564,.492,.441,.398,.360,.325,
     1  .997,.971,.852,.709,.595,.501,.435,.387,.347,.311,.279,
     2  .994,.937,.784,.634,.520,.434,.376,.333,.297,.265,.238,
     3  .985,.877,.694,.548,.442,.367,.316,.278,.247,.220,.197,
     4  .950,.773,.578,.444,.351,.289,.249,.218,.193,.172,.154,
     5  .780,.586,.415,.308,.243,.199,.171,.149,.131,.117,.105
     6  /
C
C
        KAPPA2=0.0
        IRAD=nint((100*X-25.)/5.)
c	print *,' irad = ',irad
        IF (IRAD.EQ.15)GO TO 5
        SINPHI=SIN(PHIALP)
        IF (SINPHI-.1)10,10,20
C
C  SINPHI LESS THAN .1
C
10      IF (SINPHI-.05)11,11,12
11      KAPPA2=F(1,IRAD)
        GO TO 5
12      KAPPA2=(F(2,IRAD)-F(1,IRAD))*(SINPHI-.05)/.05 + F(1,IRAD)
        GO TO 5
20      ISINPH=IFIX(SINPHI*10.)+1.
        KAPPA2=(F(ISINPH+1,IRAD)-F(ISINPH,IRAD))*(SINPHI-FLOAT(
     1  ISINPH-1)/10.)/.1 + F(ISINPH,IRAD)
5       CONTINUE
        RETURN
        END

        FUNCTION KAPPA3(X,PHIALP)
C
C   RETURNS THE GOLDSTEIN KAPPA FACTOR FOR TWO BLADED PROPELLERS
C
        REAL KAPPA3
        DIMENSION F(11,14)
        DATA F/
     1  1.00,1.00,.997,.992,.984,.975,.973,.983,1.01,1.06,1.18,
     2  1.00,1.00,.996,.990,.978,.963,.952,.949,.956,.984,1.07,
     3  1.00,.999,.996,.989,.972,.952,.930,.914,.903,.908,.950,
     4  1.00,.999,.995,.987,.966,.940,.909,.880,.849,.831,.832,
     5  1.00,.999,.993,.976,.945,.905,.863,.825,.786,.761,.751,
     6  1.00,.999,.990,.966,.923,.871,.817,.769,.724,.692,.669,
     7  1.00,.999,.988,.955,.902,.836,.771,.714,.661,.622,.588,
     8  1.00,.999,.976,.924,.856,.781,.711,.650,.597,.558,.523,
     9  1.00,.998,.964,.892,.809,.725,.650,.586,.533,.494,.457,
     1  .999,.994,.935,.843,.746,.658,.582,.520,.470,.429,.393,
     2  .998,.980,.884,.774,.670,.581,.508,.450,.404,.367,.335,
     3  .990,.948,.810,.684,.581,.496,.429,.377,.337,.304,.278,
     4  .973,.872,.693,.566,.471,.396,.341,.299,.265,.239,.218,
     5  .863,.692,.512,.406,.331,.275,.236,.203,.182,.164,.149
     6  /
C
C
        KAPPA3=0.0
        IRAD=nint((100*X-25.)/5.)
        IF (IRAD.EQ.15)GO TO 5
        SINPHI=SIN(PHIALP)
        IF (SINPHI-.1)10,10,20
C
C  SINPHI LESS THAN .1
C
10      IF (SINPHI-.05)11,11,12
11      KAPPA3=F(1,IRAD)
        GO TO 5
12      KAPPA3=(F(2,IRAD)-F(1,IRAD))*(SINPHI-.05)/.05 + F(1,IRAD)
        GO TO 5
20      ISINPH=IFIX(SINPHI*10.)+1.
        KAPPA3=(F(ISINPH+1,IRAD)-F(ISINPH,IRAD))*(SINPHI-FLOAT(
     1  ISINPH-1)/10.)/.1 + F(ISINPH,IRAD)
5       CONTINUE
        RETURN
        END


        FUNCTION KAPPA4(X,PHIALP)
C
C   RETURNS THE GOLDSTEIN KAPPA FACTOR FOR TWO BLADED PROPELLERS
C
        REAL KAPPA4
        DIMENSION F(11,14)
        DATA F/
     1  1.00,1.00,.998,.996,.991,.985,.986,.999,1.03,1.10,1.25,
     2  1.00,1.00,.998,.995,.989,.980,.975,.979,.998,1.04,1.14,
     3  1.00,1.00,.997,.995,.986,.976,.965,.959,.961,.979,1.03,
     4  1.00,1.00,.997,.994,.984,.971,.954,.939,.924,.917,.923,
     5  1.00,1.00,.996,.991,.973,.949,.921,.894,.868,.849,.842,
     6  1.00,1.00,.996,.987,.961,.927,.888,.848,.812,.781,.762,
     7  1.00,1.00,.995,.984,.950,.905,.855,.803,.756,.713,.681,
     8  1.00,1.00,.992,.965,.917,.859,.800,.743,.692,.647,.612,
     9  1.00,1.00,.989,.945,.883,.812,.745,.682,.627,.581,.543,
     1  1.00,.999,.973,.909,.830,.750,.678,.615,.562,.516,.477,
     2  .999,.993,.940,.852,.759,.674,.601,.541,.490,.447,.412,
     3  .998,.985,.882,.774,.671,.585,.517,.459,.413,.375,.345,
     4  .995,.943,.777,.651,.554,.476,.414,.369,.329,.298,.272,
     5  .945,.770,.590,.476,.396,.334,.290,.255,.228,.205,.187
     6  /
C
C
        KAPPA4=0.0
        IRAD=nint((100*X-25.)/5.)
        IF (IRAD.EQ.15)GO TO 5
        SINPHI=SIN(PHIALP)
        IF (SINPHI-.1)10,10,20
C
C  SINPHI LESS THAN .1
C
10      IF (SINPHI-.05)11,11,12
11      KAPPA4=F(1,IRAD)
        GO TO 5
12      KAPPA4=(F(2,IRAD)-F(1,IRAD))*(SINPHI-.05)/.05 + F(1,IRAD)
        GO TO 5
20      ISINPH=IFIX(SINPHI*10.)+1.
        KAPPA4=(F(ISINPH+1,IRAD)-F(ISINPH,IRAD))*(SINPHI-FLOAT(
     1  ISINPH-1)/10.)/.1 + F(ISINPH,IRAD)
5       CONTINUE
        RETURN
        END
