C*********************************************************************** 0001 C 0002 C----------------------------------------------------------------------- 0003 C PROGRAM LATDER 0004 C----------------------------------------------------------------------- 0005 C 0006 C 0007 C GIVEN VALUES OF THE AIRCRAFT GEOMETRY AND OTHER PERTINENT DATA, 0008 C THIS PROGRAM PERFORMS THE FOLLOWING: 0009 C 0010 C 1) CALCULATE NON-DIMENSIONAL STABILITY DERIVATIVES 0011 C 2) CALCULATE DIMENSIONAL STABILITY DERIVATIVES 0012 C 0013 C THE DERIVATION OF THE EQUATIONS OF MOTION ON WHICH THIS ANALYSIS IS 0014 C BASED WAS TAKEN FROM 'DYNAMICS OF THE AIRFRAME', BUREAU OF 0015 C AERONAUTICS REPORT, AE-61-4II. 0016 C 0017 C THE ANALYSIS DESCRIBED ABOVE MUST MEET THE ASSUMPTIONS IMPOSED ON 0018 C THE EQUATIONS OF MOTION WHEN THEY WERE DERIVED. 0019 C THESE ASSUMPTIONS ARE: 0020 C 0021 C 1) THE AIRFRAME IS ASSUMED TO BE A RIGID BODY. 0022 C 0023 C 2) THE EARTH IS ASSUMED TO BE FIXED IN SPACE, AND, UNLESS 0024 C SPECIFICALLY STATED OTHERWISE, THE EARTH'S ATMOSPHERE IS 0025 C ASSUMED TO BE FIXED WITH RESPECT TO THE EARTH. 0026 C 0027 C 3) THE MASS OF THE AIRPLANE IS ASSUMED TO REMAIN CONSTANT FOR 0028 C THE DURATION OF ANY PARTICULAR DYNAMIC ANALYSIS. 0029 C 0030 C 4) THE X-Z PLANE IS ASSUMED TO BE A PLANE OF SYMMETRY. 0031 C 0032 C 5) THE DISTURBANCES FROM THE STEADY FLIGHT CONDITION ARE ASSUMED 0033 C TO BE SMALL ENOUGH SO THAT THE PRODUCTS AND SQUARES OF THE 0034 C CHANGES IN VELOCITIES ARE NEGLIGIBLE IN COMPARISON WITH THE 0035 C CHANGES THEMSELVES. ALSO, THE DISTURBANCE ANGLES ARE ASSUMED 0036 C TO BE SMALL ENOUGH SO THAT THE SINES OF THESE ANGLES MAY BE 0037 C SET EQUAL TO THE ANGLES AND THE COSINES SET EQUAL TO ONE. 0038 C PRODUCTS OF THESE ANGLES ARE ALSO APPROXIMATELY ZERO AND CAN 0039 C BE NEGLECTED. AND, SINCE THE DISTURBANCES ARE SMALL, THE 0040 C CHANGE IN AIR DENSITY ENCOUNTERED BY THE AIRPLANE DURING ANY 0041 C DISTURBANCE CAN BE CONSIDERED TO BE ZERO. 0042 C 0043 C 6) DURING THE STEADY FLIGHT CONDITION, THE AIRPLANE IS ASSUMED 0044 C TO BE FLYING WITH WINGS LEVEL AND ALL COMPONENTS OF VELOCITY 0045 C ZERO EXCEPT U SUB 0. W SUB 0 = 0 BECAUSE THE STABILITY AXES 0046 C WERE CHOSEN AS THE REFERENCE AXES. 0047 C 0048 C 7) THE FLOW IS ASSUMED TO BE QUASI-STEADY. 0049 C 0050 C 0051 C 0052 C THE PERTINENT AIRPLANE CHARACTERISTICS ARE DEFINED AS FOLLOWS: 0053 C 0054 C GRAVITY (ASSUME = 32.2 FT/SEC**2 FOR THIS ALTITUDE RANGE) AND THE 0055 C GCOSGM AND GSINGM ARE THE PRODUCTS OF THE ACCELERATION DUE TO 0056 C COSINE AND SINE RESPECTIVELY OF THE INITIAL FLIGHT PATH ANGLE, 0057 C GAMMA, (USUALLY ZERO FOR LEVEL FLIGHT). 0058 C 0059 C IZZ IS THE MOMENT OF INERTIA ABOUT THE Z AXIS (FT-LBS-SEC*SEC) 0060 C 0061 C IXZ IS THE PRODUCT OF INERTIA (FT-LBS-SEC*SEC) 0062 C 0063 C IXX IS THE MOMENT OF INERTIA ABOUT THE X AXIS (FT-LBS-SEC*SEC) 0064 C 0065 C MS IS THE MASS OF THE AIRPLANE (SLUGS) 0066 C 0067 C S IS THE WING AREA OF THE AIRPLANE (SQUARE FEET) 0068 C 0069 C RHO IS THE DENSITY AT THE ALTITUDE AT WHICH THE AIRPLANE IS FLYING 0070 C 0071 C U IS THE SPEED OF THE AIRCRAFT IN FEET PER SECOND 0072 C 0073 C CH IS THE MEAN AERODYNAMIC CHORD OF THE WING (FEET) 0074 C 0075 C B IS THE WING SPAN (FEET) 0076 C 0077 C CL IS THE AIRPLANE LIFT COEFFICIENT 0078 C 0079 C SA IS THE WING SWEEP ANGLE (POSITIVE AFT, IN RADIANS) 0080 C 0081 C DIH IS THE WING DIHEDRAL ANGLE (POSITIVE UP, IN DEGREES) 0082 C 0083 C ZW IS THE DISTANCE FROM BODY CENTERLINE TO QUARTER-CHORD POINT OF 0084 C EXPOSED WING ROOT CHORD (POSITIVE FOR QUARTER-CHORD POINT BELOW 0085 C THE BODY CENTERLINE, FEET) 0086 C 0087 C FUSVOL IS THE VOLUME OF THE FUSELAGE (CUBIC FEET) 0088 C 0089 C H IS THE MAXIMUM BODY HEIGHT AT WING-BODY INTERSECTION (FEET) 0090 C 0091 C SV IS THE AREA OF THE VERTICAL TAIL(SQUARE FEET) 0092 C 0093 C BV IS THE SPAN OF THE VERTICAL TAIL(FEET) 0094 C 0095 C R1 IS THE RADIUS OF THE FUSELAGE IN THE VICINITY OF THE VERTICAL 0096 C TAIL(FEET) 0097 C 0098 C TR IS THE WING TAPER RATIO (TIP CHORD/ROOT CHORD) 0099 C 0100 C ZV IS THE DISTANCE FROM THE CENTER OF PRESSURE OF THE VERTICAL 0101 C TAIL TO THE AIRPLANE'S X-AXIS (POSITIVE FOR VERTICAL TAIL ABOVE 0102 C THE X-AXIS, FEET) 0103 C 0104 C ETAV IS THE EFFICIENCY FACTOR OF THE VERTICAL TAIL 0105 C 0106 C SBS IS THE BODY SIDE AREA OF THE FUSELAGE(SQUARE FEET) 0107 C 0108 C LF IS THE LENGTH OF THE FUSELAGE(FEET) 0109 C 0110 C LT IS THE LENGTH FROM C.G. TO CENTER OF PRESSURE OF THE TAIL, FEET 0111 C 0112 C XM IS THE DISTANCE FROM THE NOSE TO THE C.G.(FEET) 0113 C 0114 C H1 IS THE FUSELAGE HEIGHT MEASURED AT 1/4 LF FROM THE NOSE(FEET) 0115 C 0116 C H2 IS THE FUSELAGE HEIGHT MEASURED AT 3/4 LF FROM THE NOSE(FEET) 0117 C 0118 C W IS THE MAXIMUM WIDTH OF THE FUSELAGE(FEET) 0119 C 0120 C SAH IS THE SWEEP ANGLE OF THE HORIZONTAL TAIL(IN RADIANS) 0121 C 0122 C CLA2DW IS TWO-DIMENSIONAL LIFT CURVE SLOPE OF THE WING(PER RADIAN) 0123 C 0124 C BH IS THE SPAN OF THE HORIZONTAL TAIL, FEET 0125 C 0126 C SH IS THE AREA OF THE HORIZONTAL TAIL, SQUARE FEET 0127 C 0128 C TRH IS TAPER RATIO OF THE HORIZONTAL TAIL (TIP CHORD/ROOT CHORD) 0129 C 0130 C CLA2DH IS THE TWO-DIMENSIONAL LIFT CURVE SLOPE OF THE HORIZONTAL 0131 C TAIL (PER RADIAN) 0132 C 0133 C BA IS THE SPAN OF AN AILERON(ON ONE SIDE), FEET 0134 C 0135 C CA IS THE AILERON CHORD, FEET 0136 C 0137 C SR IS THE AREA OF THE RUDDER, SQUARE FEET 0138 C 0139 C ALPHA IS THE ANGLE OF ATTACK AT WHICH THE AIRPLANE IS OPERATING(IN 0140 C RADIANS) 0141 C 0142 C CDO IS THE PARASITE DRAG OF THE AIRCRAFT 0143 C 0144 C YI IS THE DISTANCE FROM THE BODY CENTERLINE TO THE INBOARD EDGE 0145 C OF THE AILERON, FEET 0146 C 0147 C THE FOLLOWING CARDS DEFINE VARIABLES USED TO APPROXIMATE FUSELAGE 0148 C VOLUME BY DIVIDING THE VOLUME INTO 4 PRISMOIDS. 0149 C 0150 C HNOSE IS THE FUSELAGE HEIGHT IN THE NOSE REGION,FEET 0151 C 0152 C WNOSE IS THE FUSELAGE WIDTH IN THE NOSE REGION,FEET 0153 C 0154 C HFCY IS THE FUSELAGE HEIGHT AT THE FRONT OF THE CANOPY, FEET 0155 C 0156 C WFCY IS THE FUSELAGE WIDTH AT THE FRONT OF THE CANOPY, FEET 0157 C 0158 C LFCY IS THE LENGTH ALONG THE BODY CENTERLINE FROM NOSE TO FRONT OF 0159 C CANOPY, FEET 0160 C 0161 C LMH IS THE LENGTH ALONG THE BODY CENTERLINE FROM NOSE TO POINT OF 0162 C MAXIMUM FUSELAGE HEIGHT, FEET 0163 C 0164 C HBCY IS THE FUSELAGE HEIGHT AT THE BACK OF THE CANOPY, FEET 0165 C 0166 C WBCY IS THE FUSELAGE WIDTH AT THE BACK OF THE CANOPY, FEET 0167 C 0168 C LBCY IS THE LENGTH ALONG THE BODY CENTERLINE FROM NOSE TO BACK OF 0169 C CANOPY, FEET 0170 C 0171 C 0172 C ALL COMPUTATIONS ARE CARRIED OUT IN DOUBLE PRECISION. 0173 C 0174 C 0175 C DATE OF CURRENT VERSION: MAY 1988 0176 C 0177 C 0178 C THIS VERSION HAS BEEN MODIFIED TO COMPILE UNDER MICROSOFT C FORTRAN77 VERSION 3.20. CHANGES NECESSARY TO PERMIT THIS C ARE GENERALLY INDICATED BY THE ABSENCE OF LINE NUMBERS. C C DATE OF MODIFICATION: JUNE 1984 C C C*********************************************************************** 0179 C IMPLICIT REAL*8(A-H,O-Z) 0182 REAL*8 MS,KGAIN,KC,KROOT,KKK,KD,KN,IXZ,IXX,IZZ,LP,NR,NP,LR,LB,NB,L 0183 1IN,NIN,NPHI,NPSI,NS,LV,NV,KY,KI,LF,LT,KB,MU,LFCY,LMH,LBCY 0184 REAL*4 AMP,W,W0 COMMON/TRANS/JPT,JEND,IVERY 0185 DIMENSION NS(5),NPHI(5),NPSI(5),DS(6),ROOTR(10),ROOTI(10),C(5),KC( 0186 14),QUF1(3),QUF2(3),RR1(2),RR2(2),RI1(2),RI2(2),CC(4),RR(3),COFFI(6 0187 1),RI(3),XR(3),XI(3),COF(3),RE(2),RIM(2),ROOT(1),COFF(2),KKK(21),KE 0188 1Y(15) 0189 C 0190 C READ IN PERTINENT AIRPLANE CHARACTERISTICS 0191 C 0192 OPEN(1,FILE='LADATA',STATUS='OLD') OPEN(3,FILE='LATDER.TXT',STATUS='NEW') WRITE(*,9000) 9000 FORMAT(1X,'invoking program LATDER') WRITE(*,9001) 9001 FORMAT(2X,'indicate how many points on the root locus you wish',/ 1,4x,'to compute') READ(*,*) NJPT IF(NJPT.LE.1) GO TO 228 WRITE(*,9002) 9002 FORMAT(5X,'indicate the number of the parameter to be varied') READ(*,*) IVERY 228 JEND=0 JPT=0 229 JPT=JPT+1 IF(JPT.GT.1) GO TO 35 0193 C----------------------------------------------------------------------- 0194 READ(1,1) CL, 0195 2 SA,DIH,ZW,H,S,SV,BV,R1,TR,ZV,ETAV,B,LF,LT,XM,H1,H2,W 0196 3,CLA2DW,SAH,BH,SH,TRH,CLA2DH,BA,CA,SR,ALPHA,CDO,RHO,YI,U,MS,GCOSGM 0197 4,GSINGM,IXX,IXZ,IZZ,HNOSE,WNOSE,HFCY,WFCY,LFCY,LMH,HBCY,WBCY,LBCY 1 FORMAT (8F10.6/8F10.6/8F10.6/8F10.6/6F10.6,F10.5,F10.6/8F10.6) 0199 C----------------------------------------------------------------------- 0200 IF(IZZ.EQ.0.0D0) IZZ=2000.0D0 ICZZ=IZZ WRITE(*,8000) WRITE(*,5050) IZZ 5050 FORMAT(7X,'IZZ=',D23.16) 8000 FORMAT(1X,'input data read successfully') 35 CONTINUE 0201 C 0202 C THE FOLLOWING CARDS PERMIT THE USER TO VARY ONE OF THE 0203 C PARAMETERS OVER SEVERAL ORDERS OF MAGNITUDE IN ORDER TO 0204 C INVESTIGATE THE EFFECT ON THE SYSTEM ROOTS. AT THE 0205 C CONCLUSION OF THE VARIATION, A PLOT OF THE ROOT LOCUS IS 0206 C PRODUCED. 0207 C IF(NJPT.LE.1.AND.JPT.GT.1) GO TO 40 0208 IF(JPT.LE.1) WW=W 0209 IF(JPT.GT.1) W=WW 0210 IF(JPT.LE.1) HH=H 0211 IF(JPT.GT.1) H=HH 0212 IF(JPT.EQ.1) GO TO 40 0213 IF(JPT.EQ.2) V=0.1D0 0214 IF(JPT.GE.3) V=1.3895D0 0215 GO TO (301,302,303,304,305,306,307,308,309,310,311,312,313, 0216 1314,315,316,317,318,319,320,321),IVERY 0217 301 IZZ=IZZ*V GO TO 40 0219 302 IXZ=IXZ*V 0220 GO TO 40 0221 303 IXX=IXX*V 0222 GO TO 40 0223 304 MS=MS*V 0224 GO TO 40 0225 305 S=S*V 0226 GO TO 40 0227 306 RHO=RHO*V 0228 GO TO 40 0229 307 U=U*V 0230 GO TO 40 0231 308 B=B*V 0232 GO TO 40 0233 309 CL=CL*V 0234 GO TO 40 0235 310 SA=SA*V 0236 GO TO 40 0237 311 DIH=DIH*V 0238 GO TO 40 0239 312 SV=SV*V 0240 GO TO 40 0241 313 BV=BV*V 0242 GO TO 40 0243 314 TR=TR*V 0244 GO TO 40 0245 315 ETAV=ETAV*V 0246 GO TO 40 0247 316 SBS=SBS*V 0248 GO TO 40 0249 317 LF=LF*V 0250 GO TO 40 0251 318 LT=LT*V 0252 GO TO 40 0253 319 CLA2DH=CLA2DH*V 0254 GO TO 40 0255 320 ALPHA=ALPHA*V 0256 GO TO 40 0257 321 CD0=CD0*V 0258 40 CONTINUE 0259 C 0260 C CALCULATE EFFECTIVE ASPECT RATIO OF VERTICAL TAIL AND LIFT CURVE 0261 C SLOPE OF VERTICAL TAIL 0262 C 0263 AE=1.55*BV*BV/SV 0264 IF (AE.GE.0.0.AND.AE.LE.6.5) GO TO 3 0265 WRITE (3,2) AE 0266 2 FORMAT (1X,'VALUE OF AE = ',G11.4,' USED TO CALCULATE CLAVT IS O 0267 1UTSIDE PREFERRED RANGE OF 0.0 TO 6.5') 0268 3 CLAVT=.0000397694*AE**5-.00069754*AE**4+.00461464*AE**3-.0158634*A 0269 1E*AE+.0401979*AE+.00037328 0270 CLAVT=CLAVT*57.2958 0271 C 0272 C CALCULATE ASPECT RATIO OF WING AND HORIZONTAL TAIL 0273 C 0274 AR=B*B/S 0275 ARH=BH*BH/SH 0276 C 0277 C CALCULATE THE MEAN AERODYNAMIC CHORD 0278 C 0279 CH=S/B 0280 C 0281 C ESTIMATE BODY SIDE AREA USING FOUR TRAPEZOIDS 0282 C 0283 SBS=(HNOSE+HFCY)*LFCY/2.+(H+HFCY)*(LMH-LFCY)/2.+(H+HBCY)*(LBCY-LMH 0284 1)/2.+(HBCY+2.*R1)*(LF-LBCY)/2. 0285 C 0286 C ESTIMATE FUSELAGE VOLUME USING FOUR PRISMOIDS 0287 C 0288 V1=LFCY*(2.0*(HNOSE*WNOSE+HFCY*WFCY)+HFCY*WNOSE+HNOSE*WFCY) 0289 V2=(LMH-LFCY)*(2.0*(HFCY*WFCY+H*W)+H*WFCY+HFCY*W) 0290 V3=(LBCY-LMH)*(2.0*(HBCY*WBCY+H*W)+H*WBCY+W*HBCY) 0291 V4=(LF-LBCY)*(2.0*(HBCY*WBCY+(2.0*R1)*(2.0*R1))+HBCY*2.0*R1+WBCY*2 0292 1.0*R1) 0293 FUSVOL=(V1+V2+V3+V4)/6.0 0294 C 0295 C NEXT THREE CARDS DETERMINE WHETHER THE AIRPLANE IS LOW WING(WP=3.0 0296 C MID WING(WP=2.0), OR HIGH WING(WP=1.0) 0297 C 0298 R=ZW/H 0299 WP=2.0 0300 IF (R.GE..25)WP=3.0 0301 IF (R.LT.-.25)WP=1.0 0302 IF(IZZ.EQ.0.0D0) IZZ=ICZZ C 0303 C WRITE PERTINENT AIRPLANE CHARACTERISTICS AS DEFINED IN THE INITIAL 0304 C PART OF THIS PROGRAM 0305 C 0306 WRITE (3,4) RHO,S,MS,GCOSGM,U,CH,B,GSINGM,IXX,IXZ,IZZ,CL,SA,DIH,ZW 1,FUSVOL,H,SV,BV,R1,TR,ZV,ETAV,SBS,LF,LT,XM,H1,H2,W,SAH,CLA2DW,BH,S 0308 1H,TRH,CLA2DH,BA,CA,SR,ALPHA,CDO,YI,HNOSE,WNOSE 0309 4 FORMAT (1X,121('*'),/,1X,'*',119X,'*',/,1X,'*',40X,'PERTINENT AIRP 0310 1LANE CHARACTERISTICS',45X,'*',/,1X,'*',40X,34('-'),45X,'*',/,1X,'* 0311 1',119X,'*',/,1X,'*',17X,'RHO =',F10.6,3X,'WING AREA =',F10.6,3X,'M 0312 1ASS =',F10.6,3X,'G*COS(GAMMA) =',F10.6,17X,'*',/,1X,'*',17X,'U = 0313 1',F10.4,3X,'CHORD =',F10.6,3X,'SPAN =',F10.4,3X,'G*SIN(GAMMA) 0314 1=',F10.6,17X,'*'/,1X,'*',17X,'IXX =',F10.4,3X,'IXZ =',F9.4,4 0315 1X,'IZZ =',F9.4,4X,'CL',11X,'=',F10.6,17X,'*',/,1X,'*',17X,'SA =' 0316 1,F10.6,3X,'DIH',7X,'=',F10.6,3X,'ZW =',F10.6,3X,'FUSVOL',7X,'=', 0317 1F10.6,17X,'*',/,1X,'*',17X,'H =',F10.6,3X,'SV',8X,'=',F10.6,3X,' 0318 1BV =',F10.6,3X,'R1',11X,'=',F10.6,17X,'*',/,1X,'*',17X,'TR =',F 0319 110.6,3X,'ZV',8X,'=',F10.6,3X,'ETAV =',F10.6,3X,'SBS',10X,'=',F10.6 0320 1,17X,'*',/,1X,'*',17X,'LF =',F10.6,3X,'LT',8X,'=',F10.6,3X,'XM 0321 1=',F10.6,3X,'H1',11X,'=',F10.6,17X,'*',/,1X,'*',17X,'H2 =',F10.6, 0322 13X,'W',9X,'=',F10.6,3X,'SAH =',F10.6,3X,'CLA2DW',7X,'=',F10.6,17X 0323 1,'*',/,1X,'*',17X,'BH =',F10.6,3X,'SH',8X,'=',F10.6,3X,'TRH =',F 0324 110.6,3X,'CLA2DH',7X,'=',F10.6,17X,'*',/,1X,'*',17X,'BA =',F10.6,3 0325 1X,'CA',8X,'=',F10.6,3X,'SR =',F10.6,3X,'ALPHA',8X,'=',F10.6,17X, 0326 1'*',/,1X,'*',17X,'CDO =',F10.6,3X,'YI',8X,'=',F10.6,3X,'HNOSE=',F1 0327 10.6,3X,'WNOSE',8X,'=',F10.6,17X,'*') 0328 WRITE (3,5) HFCY,WFCY,LFCY,LMH,HBCY,WBCY,LBCY 0329 5 FORMAT (1X,'*',17X,'HFCY=',F10.6,3X,'WFCY',6X,'=',F10.6,3X,'LFCY = 0330 1',F10.6,3X,'LMH',10X,'=',F10.6,17X,'*',/,1X,'*',17X,'HBCY=',F10.6, 0331 13X,'WBCY',6X,'=',F10.6,3X,'LBCY =',F10.6,44X,'*',/,1X,'*',119X,'*' 0332 1,/,1X,121('*'),/) 0333 C 0334 C THE FOLLOWING READ STATEMENT ENABLES ONE TO READ IN A VALUE FOR A 0335 C PARTICULAR STABILITY DERIVATIVE OR GROUP OF DERIVATIVES, THEREBY 0336 C OVER RIDING THE CALCULATED VALUE, BY SIMPLY PUTTING A ONE (1.0) IN 0337 C COLUMN 5 FOR CYB, 10 FOR CLB, 15 FOR CNB, 20 FOR CLP, 25 FOR CYP, 0338 C 30 FOR CNP, 35 FOR CYR, 40 FOR CLR, 45 FOR CNR, 50 FOR CYDA, 55 0339 C FOR CLDA, 60 FOR CNDA, 65 FOR CYDR, 70 FOR CLDR, OR 75 FOR CNDR, 0340 C AND THEN PLACING THE VALUES TO BE READ, IN ORDER, (ONE PER CARD) 0341 C BEHIND THIS KEY CARD. 0342 C 0343 IF(JPT.GT.1) GO TO 36 0344 C----------------------------------------------------------------------- 0345 READ (1,6) (KEY(JI),JI=1,15) 0346 C----------------------------------------------------------------------- 0347 WRITE(*,8001) WRITE(*,6) (KEY(JI),JI=1,15) 8001 FORMAT(3X,'stability derivative key selection values read in') 6 FORMAT (15I5) 0348 36 CONTINUE 0349 C 0350 C THE FOLLOWING CARDS CALL PROGRAM SECTIONS WHICH COMPUTE 0351 C VALUES FOR THE NON-DIMENSIONAL STABILITY DERIVATIVES 0352 C 0353 GO TO 1001 0354 901 IF (KEY(1).EQ.0) GO TO 8 0355 READ (1,7) CYB 0356 7 FORMAT (F10.6) 0357 8 GO TO 1002 0358 902 IF (KEY(2).EQ.0) GO TO 9 0359 READ (1,7) CLB 0360 9 GO TO 1003 0361 903 IF (KEY(3).EQ.0) GO TO 10 0362 READ (1,7) CNB 0363 10 GO TO 1004 0364 904 IF (KEY(4).EQ.0) GO TO 11 0365 READ (1,7) CLP 0366 11 GO TO 1005 0367 905 IF (KEY(5).EQ.0) GO TO 12 0368 READ (1,7) CYP 0369 12 GO TO 1006 0370 906 IF (KEY(6).EQ.0) GO TO 13 0371 READ (1,7) CNP 0372 13 GO TO 1007 0373 907 IF (KEY(7).EQ.0) GO TO 14 0374 READ (1,7) CYR 0375 14 GO TO 1008 0376 908 IF (KEY(8).EQ.0) GO TO 15 0377 READ (1,7) CLR 0378 15 GO TO 1009 0379 909 IF (KEY(9).EQ.0) GO TO 16 0380 READ (1,7) CNR 0381 16 GO TO 1010 0382 910 IF (KEY(10).EQ.0) GO TO 17 0383 READ (1,7) CYDA 0384 17 GO TO 1011 0385 911 IF (KEY(11).EQ.0) GO TO 18 0386 READ (1,7) CLDA 0387 18 GO TO 1012 0388 912 IF (KEY(12).EQ.0) GO TO 19 0389 READ (1,7) CNDA 0390 19 GO TO 1013 0391 913 IF (KEY(13).EQ.0) GO TO 20 0392 READ (1,7) CYDR 0393 20 GO TO 1014 0394 914 IF (KEY(14).EQ.0) GO TO 21 0395 READ (1,7) CLDR 0396 21 GO TO 1015 0397 915 IF (KEY(15).EQ.0) GO TO 22 0398 READ (1,7) CNDR 0399 22 WRITE (3,23) 0400 23 FORMAT (1X,121('*'),/,1X,'*',119X,'*',/,1X,'*',44X,'LATERAL STABIL 0401 1ITY DERIVATIVES',46X,'*',/,1X,'*',44X,29('-'),46X,'*',/,1X,'*',119 0402 1X,'*',/,1X,'*',119X,'*') 0403 WRITE (3,24) CYB,CLB,CNB,CYP,CLP,CNP,CYR,CLR,CNR,CYDA,CLDA,CNDA,CY 0404 1DR,CLDR,CNDR 0405 24 FORMAT (1X,'* CYB =',F10.6,4X,'CLB =',F10.6,4X,'CNB =',F10.6,4 0406 1X,'CYP =',F10.6,4X,'CLP =',F10.6,4X,'CNP =',F10.6,' *',/,1X,'* 0407 1 CYR =',F10.6,4X,'CLR =',F10.6,4X,'CNR =',F10.6,4X'CYDA =',F9.6 0408 1,4X,'CLDA =',F10.6,4X,'CNDA =',F10.6,' *',/,1X,'* CYDR =',F10.6, 0409 14X,'CLDR =',F10.6,4X,'CNDR =',F9.6,59X,' *',/,1X,'*',119X,'*',/, 0410 11X,121('*'),/) 0411 C 0412 WRITE(*,8002) 8002 FORMAT(4X,'derivative values now computed or read in') IF(JPT.GT.1) GO TO 37 0413 C THE READ STATEMENT BELOW IS THE INPUT FOR THE VARIABLES EXPLAINED 0414 C BELOW. 0415 C---------------------------------------------------------------------- 0416 25 READ (1,26) K,NUMER,ACC 0417 C---------------------------------------------------------------------- 0418 WRITE(*,26) K,NUMER,ACC 26 FORMAT (2I3,E14.7) 0419 37 CONTINUE 0420 C 0421 IF (K.EQ.1) GO TO 27 0422 CYIN=CYDA 0423 CLIN=CLDA 0424 CNIN=CNDA 0425 GO TO 28 0426 27 CYIN=CYDR 0427 CLIN=CLDR 0428 CNIN=CNDR 0429 C 0430 C CYIN, CLIN, AND CNIN ARE THE STABILITY DERIVATIVES DUE TO CONTROL 0431 C SURFACE DEFLECTIONS. CYIN IS EITHER THE PARTIAL OF CY WITH 0432 C RESPECT TO RUDDER DEFLECTION OR AILERON DEFLECTION. THE VALUE OF K 0433 C DETERMINES WHETHER IT IS CONCERNED WITH RUDDER OR AILERON. 0434 C 0435 C IF K IS GIVEN THE VALUE 1 THEN CYIN, CLIN, AND CNIN ARE PARTIAL 0436 C DERIVATIVES WITH RESPECT TO RUDDER DEFLECTION. IF K=2 THEN THE 0437 C PARTIALS ARE TAKEN WITH RESPECT TO AILERON DEFLECTION. 0438 C 0439 C ACC IS THE ACCURACY USED THROUGHOUT THE PROGRAM TO COMPARE TWO 0440 C NUMBERS TO SEE IF THEY ARE SUFFICIENTLY CLOSE. 0441 C 0442 28 IF (K.NE.1) GO TO 30 0443 WRITE (3,29) 0444 29 FORMAT (15X,98('*'),/,15X,'*',96X,'*',/,15X,'* ',32X,'RESPONSE TO 0445 1 RUDDER DEFLECTION',31X,' *') 0446 GO TO 32 0447 30 WRITE (3,31) 0448 31 FORMAT (15X,98('*'),/,15X,'*',96X,'*',/,15X,'* ',31X,'RESPONSE TO 0449 1 AILERON DEFLECTION',31X,' *') 0450 32 WRITE (3,33) CYIN,CLIN,CNIN,K,ACC 0451 33 FORMAT (15X,'*',96X,'*',/,15X,'*',96X,'*',/15X,'* ','CYIN =',F10. 0452 16,4X,'CLIN =',F10.6,4X,'CNIN =',F10.6,4X,'K =',I2,4X,'ACC =',F14.8 0453 1,6X,'*',/,15X,'*',96X,'*',/,15X,98('*')) 0454 C 0455 C CALCULATION OF DIMENSIONAL STABILITY DERIVATIVES 0456 C 0457 A=RHO*U*S 0458 F=RHO*U*S*B 0459 H=RHO*U*S*B*B 0460 D=RHO*U*U*S 0461 E=RHO*U*U*S*B 0462 C 0463 C A, F, H, D, E ARE JUST CONSTANTS USED TO CALCULATE THE DIMENSIONAL 0464 C STABILITY DERIVATIVES. 0465 C 0466 IF(IZZ.EQ.0.0D0) IZZ=ICZZ WRITE(*,6200) MS,IXX,IZZ 6200 FORMAT(2X,'MS=',F15.8,1X,'IXX=',F15.8,1X,'IZZ=',F15.8) YV=A*CYB/(2.0D0*MS) 0467 LV=F*CLB/(2.0D0*IXX) 0468 LB=U*LV 0469 NV=F*CNB/(2.0D0*IZZ) 0470 NB=U*NV 0471 YP=F*CYP/(4.0D0*MS) 0472 LP=H*CLP/(4.0D0*IXX) 0473 NP=H*CNP/(4.0D0*IZZ) 0474 YR=F*CYR/(4.0D0*MS) 0475 LR=H*CLR/(4.0D0*IXX) 0476 NR=H*CNR/(4.0D0*IZZ) 0477 YIN=D*CYIN/(2.0D0*MS) 0478 LIN=E*CLIN/(2.0D0*IXX) 0479 NIN=E*CNIN/(2.0D0*IZZ) 0480 WRITE(*,8005) 8005 FORMAT(5X,'dimensional stability derivatives now calculated') C 0481 C THE FOLLOWING WRITE STATEMENT GIVES THE DIMENSIONAL STABILITY 0482 C DERIVATIVES. 0483 C 0484 WRITE (3,34) YV,LB,NB,YP,LP,NP,YR,LR,NR,YIN,LIN,NIN 0485 34 FORMAT (/2X,123('*'),/,2X,'*',121X,'*',/,2X,'*',44X,'DIMENSIONAL S 0486 1TABILITY DERIVATIVES',44X,'*',/,2X,'*',121X,'*',/,2X,'*',2X,' YV 0487 1=',F10.5,4X,' LB =',F10.5,4X,' NB =',F10.5,4X,' YP =',F10.5,4X, 0488 1' LP =',F10.5,4X,' NP =',F10.5,3X,'*',/,2X,'*',3X,' YR =',F10.5, 0489 14X,' LR =',F10.5,4X,' NR =',F10.5,4X,' YIN =',F10.5,4X,' LIN =', 0490 1F10.5,4X,' NIN =',F10.5,3X,'*',/,2X,'*',121X,'*',/,2X,123('*')) 0491 C 0492 IF(JPT.GT.1) GO TO 38 0493 C---------------------------------------------------------------------- 0494 READ(1,899) AMP,W,W0,IF,ISET,LK,LBB 0495 C---------------------------------------------------------------------- 0496 IF=5 WRITE(*,899) AMP,W,W0,IF,ISET,LK,LBB 899 FORMAT(3F10.5,4I10) 0497 38 CONTINUE 0498 C 0499 C IF LBB IS 6 OR GREATER BRANCH TO STATIC STABILITY AND 0500 C CONTROL COMPUTATION, FORCEL 0501 C 0502 IF(LBB.GE.6) CALL FORCEL(ETAV,SV,BV,U,RHO,MS,B,CH,BA,SR, 0503 1CA,YI,CNB,CNDA,CLDR,CLDA,CNDR,CLB,IXX,CLP,S,LBB) 0504 C 0505 IF(LBB.EQ.9) STOP 0506 C 0507 CALL LATPRG(YV,U,GCOSGM,GSINGM,YP,YR,LB,LP,LR, 0508 1IXX,IZZ,IXZ,NB,NP,NR,YIN,LIN,NIN,IF,AMP,ISET,W,W0,LK,LBB) 0509 C 0510 C*********************************************************************** 0511 C 0512 C LATDER TERMINATES AT THIS POINT. 0513 C 0514 C*********************************************************************** 0515 C 0516 IF(JPT.EQ.(NJPT-1)) JEND=1 IF(JPT.LE.(NJPT-1)) GO TO 229 STOP 0517 C 0518 C*********************************************************************** 0519 C 0520 C THE FOLLOWING SECTIONS CALCULATE THE INDIVIDUAL STABILITY 0521 C DERIVATIVE VALUES. 0522 C 0523 C*********************************************************************** 0524 C 0525 1001 CYBWS=CL*CL*6.0*DTAN(SA)*DSIN(SA)/(3.14159*AR*(AR+4.0*DCOS(SA))) 0526 CYBWD=-.0001*DABS(DIH)*57.2958 0527 CLAFUS=.1 0528 BRA=FUSVOL**.66666666666666 0529 IF (ZW.LT.0.0) GO TO 401 0530 KI=1.0+.485*(ZW*2.0/H) 0531 GO TO 402 0532 401 KI=1.0-.85555*(ZW*2.0/H) 0533 402 CYBFUS=-KI*CLAFUS*BRA/S 0534 SWDPR=.724+3.06*SV/(S*(1.0+DCOS(SA)))+.4*ZW/H+.009*AR 0535 PAR=BV/(2.0*R1) 0536 KY=1.0+.166*(PAR-3.5) 0537 IF (PAR.LE.2.0)KY=.75 0538 IF (PAR.GE.3.5)KY=1.0 0539 CYBT=-KY*CLAVT*SWDPR*SV/S 0540 CYB=CYBWS+CYBWD+CYBFUS+CYBT 0541 WRITE(*,7000) CYB 7000 FORMAT(1X,'CYB=',F15.8) GO TO 901 0542 C 0543 1002 IF (AR.GE.1.5.AND.AR.LE.10.) GO TO 412 0544 WRITE (3,411) AR 0545 411 FORMAT (1X,'VALUE OF AR = ',G11.4,2X,'USED TO CALCULATE CLB IS OUT 0546 1SIDE PREFERRED RANGE OF 1.5 TO 10.0') 0547 412 CLBOD=((-.20833-DSQRT(.0434-4.6296*(1.0052-AR)))/2.3148148+.7*(1.0 0548 1-TR)*(1.0-TR)*(AR-1.5)/6.5)*.0001 0549 CLBWD=CLBOD*DIH*57.2958 0550 CLBWWD=CL*(-1.25*(.71*TR+.29)/(AR*TR)+.05) 0551 CLBW=CLBWD+CLBWWD 0552 CLBVT=-CLAVT*SV*ZV*ETAV/(S*B) 0553 IF (WP.NE.1.0) GO TO 413 0554 DCLB=-.00044 0555 GO TO 415 0556 413 IF (WP.NE.2.0) GO TO 414 0557 DCLB=0.0 0558 GO TO 415 0559 414 DCLB=.00064 0560 415 CLB=CLBW+CLBVT+DCLB 0561 WRITE(*,7001) CLB 7001 FORMAT(1X,'CLB=',F15.8) GO TO 902 0562 C 0563 1003 A=XM/LF 0564 IF (A.GE..1.AND.A.LE..8) GO TO 422 0565 WRITE (3,421) A 0566 421 FORMAT (1X,'VALUE OF XM/LF = ',G11.4,' USED TO CALCULATE CNB IS 0567 1OUTSIDE PERFERRED RANGE OF 0.1 TO 0.8') 0568 422 X=(A+.2)*10.0 0569 XP=X-3.0 0570 R=LF*LF/SBS 0571 IF (R.GE.2.5.AND.R.LE.20.) GO TO 424 0572 WRITE (3,423) R 0573 423 FORMAT (1X,'VALUE OF LF*LF/SBS = ',G11.4,' USED TO CALCULATE CNB 0574 1 IS OUTSIDE PREFERRED RANGE OF 2.5 TO 20.0') 0575 424 YO=-.0015609*R**3+.0675772*R**2-.999884*R+1.20213 0576 YP=.314*XP+YO 0577 Y=12.0+YP 0578 YP=Y-8.0 0579 AA=DSQRT(H1/H2) 0580 IF (AA.GE..8.AND.AA.LE.1.65) GO TO 426 0581 WRITE (3,425) AA 0582 425 FORMAT (1X,'VALUE OF SQRT(H1/H2) = ',G11.4,' USED TO CALCULATE CN 0583 1B IS OUTSIDE PREFERRED RANGE OF 0.8 TO 1.65') 0584 426 XP=AA*YP 0585 AB=H/W 0586 IF (AB.GE..5.AND.AB.LE.2.0) GO TO 428 0587 WRITE (3,427) AB 0588 427 FORMAT (1X,'VALUE OF H/W = ',G11.4,' USED TO CALCULATE CNB IS OU 0589 1TSIDE PREFERRED RANGE OF 0.5 TO 2.0') 0590 428 YP=-((H/W)**(.48*W/H))*XP 0591 Y1=YP 0592 RHO1=RHO*10000. 0593 MU=.34777*RHO1+29.23345 0594 RN=RHO1*U*LF/(MU*.0001) 0595 AC=RN*.000001 0596 IF (AC.GE.0.0.AND.AC.LE.80.0) GO TO 430 0597 WRITE (3,429) AC 0598 429 FORMAT (1X,'VALUE OF RN*10-6 = ',G11.4,' USED TO CALCULATE CNB I 0599 1S OUTSIDE PREFERRED RANGE OF 0.0 TO 80.0') 0600 430 YP=AC/20.0 0601 X=.00833327*YP**3-.113571*YP*YP+1.02095*YP+.446857 0602 B1=X/2.0 0603 X=-2.0*(Y1-B1) 0604 KB=(X-1.0)*.0005 0605 CNB=-KB*SBS*LF*57.298/(S*B)-CYBT*LT/B 0606 WRITE(*,7002) CNB 7002 FORMAT(1X,'CNB=',F15.8) GO TO 903 0607 C 0608 1004 IF (SA.GE.-.524.AND.SA.LE..524) GO TO 432 0609 WRITE (3,431) SA 0610 431 FORMAT (1X,'VALUE OF SA = ',G11.4,' USED TO CALCULATE CYP IS OUT 0611 1SIDE PREFERRED RANGE OF -.524 TO .524 RADIANS') 0612 432 IF (AR.GE.5.) GO TO 434 0613 WRITE (3,433) AR 0614 433 FORMAT (1X,'VALUE OF AR = ',G11.4,' USED TO CALCULATE CYP IS LESS 0615 1 THAN THE SUGGESTED MINIMUM OF 5.0;',/,1X,'HOWEVER, DUE TO THE NEG 0616 1LIGIBLE EFFECT OF THIS DERIVATIVE, ONE MAY IGNORE THE EFFECT OF SM 0617 1ALLER AR ON CYP') 0618 434 X=-(SA*57.3-30.0)*7.0/60.0 0619 Y=X+3.0*(1.0-TR**1.58495) 0620 CYPOCL=(4.0-Y)/10.0 0621 IF (DIH.GE.-10..AND.DIH.LE.10.) GO TO 436 0622 WRITE (3,435) DIH 0623 435 FORMAT (1X,'VALUE OF DIH = ',G11.4,' USED TO CALCULATE CYP IS OU 0624 1TSIDE PREFERRED RANGE OF -10.0 TO 10.0 DEGREES') 0625 436 P=(DIH-10.)/17.5+.5 0626 CYP=CYPOCL*CL+P*CLP 0627 WRITE(*,7003) CYP 7003 FORMAT(1X,'CYP=',F15.8) GO TO 904 0628 C 0629 1005 IF (AR.GE.0.0.AND.AR.LE.10.) GO TO 442 0630 WRITE (3,441) AR 0631 441 FORMAT (1X,'VALUE OF AR = ',G11.4,' USED TO CALCULATE CLP IS OUTS 0632 1IDE PREFERRED RANGE OF 0.0 TO 10.0') 0633 442 WRITE(*,6000) AR 6000 FORMAT(5X,'AR=',F15.8) CLPAO=(.00000929805*AR**5-.00031818*AR**4+.00379426*AR**3-.0157796 0634 1*AR*AR-.0499261*AR-.0500059)+.209*(1.0-TR)**2*(AR-1.0)/9.0 0635 WRITE(*,6001) CLPAO,CLA2DW 6001 FORMAT(5X,'CLPAO=',F15.8,1X,'CLA2DW=',F15.8) IF(DABS(SA).GT.0.524D0) SA=0.0D0 CLPAOW=CLPAO*(AR+4.0D0*DCOS(SA))/(2.0D0*3.14159D0*AR/CLA2DW+4.0D0 1*DCOS(SA)) WRITE(*,6002) CLPAOW 6002 FORMAT(5X,'CLPAOW=',F15.8) CLPW=CLPAOW-(CL*CL/(8.0*3.14159*AR*DCOS(SA)**2))*(1.0+2.0*DSIN(SA) 0638 1**2*(AR+2.0*DCOS(SA))/(AR+4.0*DCOS(SA)))-CDO/8.0 0639 WRITE(*,6003) CLPW 6003 FORMAT(5X,'CLPW=',F15.8) CLPAOH=(.00000929805*ARH**5-.00031818*ARH**4+.00379426*ARH**3-.015 0640 17796*ARH**2-.0499261*ARH-.0500059)+.209*(1.0-TRH)**2*(ARH-1.0)/9.0 0641 WRITE(*,6004) CLPAOH 6004 FORMAT(5X,'CLPAOH=',F15.8) WRITE(*,6104) S,B,ARH,CLA2DH 6104 FORMAT(5X,'S=',F15.8,1X,'B=',F15.8,1X,'ARH=',F10.3,1X,'CLA2DH=', 1F10.3) CLPH=((0.5D0*SH*BH*BH)/(S*B*B))*(CLPAOH*(ARH+4.0D0*DCOS(SAH)) 1/((2.0D0*3.14159D0*ARH/CLA2DH)+4.0D0*DCOS(SAH))) CLPV=2.0*ZV*ZV*CYBT/(B*B) 0644 CLP=CLPW+CLPH+CLPV 0645 WRITE(*,7005) CLP 7005 FORMAT(1X,'CLP=',F15.8) GO TO 905 0646 C 0647 1006 IF (AR.LE.16.0) GO TO 452 0648 WRITE (3,451) AR 0649 451 FORMAT (1X,'VALUE OF AR = ',G11.4,' USED TO CALCULATE CNP IS OUT 0650 1SIDE PREFERRED RANGE OF 0.0 TO 16.0') 0651 452 CNPOCL=-.0000140468*AR**3+.000649352*AR*AR-.0122523*AR+.00192856 0652 CNPW=(CL*(AR+4.0)/(AR+4.0*DCOS(SA)))*(1.0+6.0*(1.0+DCOS(SA)/AR)*DT 0653 1AN(SA)*DTAN(SA)/12.0)*CNPOCL 0654 P=ZV*2.0/B 0655 IF (P.GE.0.0.AND.P.LE..6) GO TO 454 0656 WRITE (3,453) P 0657 453 FORMAT (1X,'VALUE OF ZV/(B/2) = ',G11.4,' USED TO CALCULATE CNP 0658 1IS OUTSIDE PREFERRED RANGE OF 0.0 TO 0.6') 0659 454 SIG1=22360.7*P**9-69947.1*P**8+93669.5*P**7-70229.7*P**6+32373.8*P 0660 1**5-9487.01*P**4+1764.41*P**3-200.25*P*P+11.9843*P-.00142786 0661 IF (P.GE..157) GO TO 455 0662 SIG2=216915.*P**9-610122.*P**8+724643.*P**7-473060.*P**6+185140.*P 0663 1**5-44463.6*P**4+6444.43*P**3-531.134*P**2+21.1205*P+.0395735 0664 GO TO 456 0665 455 SIG2=.533183*P*P-.78216*P+.35117 0666 456 SIGMA1=SIG2+TR*(SIG1-SIG2) 0667 DELHOB=(ZV-(ZV*DCOS(ALPHA)-LT*DSIN(ALPHA)))/B 0668 SIGMA2=9.3*DELHOB*DELHOB*3.0/B 0669 CNPV=(1.00*CLAVT*SV*(ZV*DSIN(ALPHA)+LT*DCOS(ALPHA))/(S*B))*(2.0*(Z 0670 1V*DCOS(ALPHA)-LT*DSIN(ALPHA))/B-(SIGMA1+SIGMA2)) 0671 CNP=CNPW+CNPV 0672 WRITE(*,7006) CNP 7006 FORMAT(1X,'CNP=',F15.8) GO TO 906 0673 C 0674 1007 CYRW=.143*CL-.05 0675 CYRT=-2.0*LT*CYBT/B 0676 CYR=CYRW+CYRT 0677 GO TO 907 0678 C 0679 1008 IF (AR.GE.2.0.AND.AR.LE.12.) GO TO 462 0680 WRITE (3,461) AR 0681 461 FORMAT (1X,'VALUE OF AR = ',G11.4,' USED TO CALCULATE CLR IS OUT 0682 1SIDE PREFERRED RANGE OF 2.0 TO 12.0') 0683 462 YP=-.0000247674*AR**4+.00194154*AR**3-.0468052*AR*AR+.456404*AR-.0 0684 1980299+2.15*(AR+32.4)*(TR**(1.12*(1.0-TR)))/44.4 0685 Y=YP+2.0 0686 IF (SA.GE.0.0.AND.SA.LE..524) GO TO 464 0687 WRITE (3,463) SA 0688 463 FORMAT (1X,'VALUE OF SA = ',G11.4,' USED TO CALCULATE CLR IS OUT 0689 1SIDE PREFERRED RANGE OF 0.0 TO 0.524 RADIANS') 0690 464 X=100.0*Y/(100.0-SA*57.3) 0691 CLROCL=X/20.0 0692 CLRW=CLROCL*CL 0693 CLRT=-2.0*LT*ZV*CYBT/(B*B) 0694 CLR=CLRW+CLRT 0695 WRITE(*,7008) CLR 7008 FORMAT(1X,'CLR=',F15.8) GO TO 908 0696 C 0697 1009 CNR=2.0*LT*LT*CYBT/(B*B)-(.33*(1.0+3.0*TR)*CDO/(2.0+2.0*TR)-.02*(1 0698 1.0-(AR-6.0)/13.0-(1.0-TR)/2.5)*CL*CL) 0699 GO TO 909 0700 C 0701 1010 CYDA=0.0 0702 GO TO 910 0703 C 0704 1011 IF (AR.GE.4.0.AND.AR.LE.12.0) GO TO 472 0705 WRITE (3,471) AR 0706 471 FORMAT (1X,'VALUE OF AR = ',G11.4,' USED TO CALCULATE CLDA IS O 0707 1UTSIDE PREFERRED RANGE OF 4.0 TO 12.0') 0708 472 K=0 0709 P=YI*2.0/B 0710 473 CLDABT=51.1702*P**8-238.294*P**7+451.162*P**6-444.985*P**5+241.797 0711 1*P**4-70.2549*P**3+10.5819*P*P-.408045*P+.0000119225 0712 CLDAAT=2.0833*P**5-6.1553*P**4+5.36932*P**3-.567235*P*P+.170341*P- 0713 1.000227237 0714 IF (AR.GT.6.0.AND.AR.LT.10.0) GO TO 474 0715 IF (AR.LE.6.0)CLDAOT=CLDABT 0716 IF (AR.GE.10.0)CLDAOT=CLDAAT 0717 GO TO 475 0718 474 CLDAOT=CLDABT+(CLDAAT-CLDABT)*(AR-6.0)/4.0 0719 475 IF (K.EQ.1) GO TO 476 0720 CLDALT=CLDAOT 0721 P=(YI+BA)*2.0/B 0722 K=K+1 0723 GO TO 473 0724 476 CLDAOT=CLDAOT-CLDALT 0725 R=CA/CH 0726 IF (R.GE.0.0.AND.R.LE..4) GO TO 478 0727 WRITE (3,477) R 0728 477 FORMAT (1X,'VALUE OF CA/CH = ',G11.4,' USED TO CALCULATE CLDA IS 0729 1 OUTSIDE PREFERRED RANGE OF 0.0 TO 0.4') 0730 478 T=-65.1062*R**4+50.8227*R**3-15.7949*R*R+3.53383*R+.000043467 0731 CLDA=CLDAOT*T 0732 WRITE(*,7011) CLDA 7011 FORMAT(1X,'CLDA=',F15.8) GO TO 911 0733 C 0734 1012 R=YI*2.0/B 0735 IF (AR.GE.3.0.AND.AR.LE.8.0) GO TO 482 0736 WRITE (3,481) AR 0737 481 FORMAT (1X,'VALUE OF AR = ',G11.4,' USED TO CALCULATE CNDA IS OU 0738 1TSIDE PREFERRED RANGE OF 3.0 TO 8.0') 0739 482 IF (AR.LT.3.0) GO TO 483 0740 IF (AR.GE.3.0.AND.AR.LE.4.0) GO TO 484 0741 IF (AR.GT.4.0.AND.AR.LT.6.0) GO TO 485 0742 IF (AR.GT.4.0.AND.AR.LE.8.0) GO TO 486 0743 EF1=-.110526*R*R+.013893*R-.146355 0744 EF2=-.191365*R**3+.147284*R*R-.00039296*R-.119884 0745 GO TO 487 0746 483 EF1=-.36 0747 EF2=-.0889976*R**3-.00666109*R*R+.0419053*R-.284843 0748 GO TO 487 0749 484 AF1=-.36 0750 BF1=.0538037*R**3-.133855*R*R+.00627854*R-.262321 0751 EF1=AF1+(AR-3.0)*(BF1-AF1) 0752 AF2=-.0889976*R**3-.00666109*R*R+.0419053*R-.28483 0753 BF2=.0549387*R**3-.164298*R*R+.101864*R-.234861 0754 EF2=AF2+(AR-3.0)*(BF2-AF2) 0755 GO TO 487 0756 485 AF1=.0538037*R**3-.133855*R*R+.00627854*R-.262321 0757 BF1=-.0629823*R**3-.0173375*R*R-.0333427*R-.180026 0758 EF1=AF1+((AR-4.0)/2.0)*(BF1-AF1) 0759 AF2=.0549387*R**3-.164298*R*R+.101864*R-.234861 0760 BF2=-.108911*R**3+.0291149*R*R+.0440961*R-.160312 0761 EF2=AF2+((AR-4.0)/2.0)*(BF2-AF2) 0762 GO TO 487 0763 486 AF1=-.0629823*R**3-.0173375*R*R-.0333427*R-.180026 0764 BF1=-.110526*R*R+.013893*R-.146355 0765 EF1=AF1+((AR-6.0)/2.0)*(BF1-AF1) 0766 AF2=-.108911*R**3+.0291149*R*R+.0440961*R-.160312 0767 BF2=-.191365*R**3+.147284*R*R-.00039296*R-.119884 0768 EF2=AF2+((AR-6.0)/2.0)*(BF2-AF2) 0769 487 EFF=EF2+((TR-.25)/.75)*(EF1-EF2) 0770 CNDA=2.0*EFF*CL*CLDA 0771 WRITE(*,7012) CNDA 7012 FORMAT(1X,'CNDA=',F15.8) GO TO 912 0772 C 0773 1013 R=SR/SV 0774 IF (R.GE.0.0.AND.R.LE..7) GO TO 492 0775 WRITE (3,491) R 0776 491 FORMAT (1X,'VALUE OF SR/SV = ',G11.4,' USED TO CALCULATE CYDR IS 0777 1 OUTSIDE PREFERRED RANGE OF 0.0 TO 0.7') 0778 492 T=21.7949*R**5-46.4744*R**4+36.9347*R**3-14.259*R*R+3.70551*R-.000 0779 1057815 0780 CYDR=CLAVT*T*SV/S 0781 WRITE(*,7013) CYDR 7013 FORMAT(1X,'CYDR=',F15.8) GO TO 913 0782 C 0783 1014 R=SR/SV 0784 IF (R.GE.0.0.AND.R.LE..7) GO TO 502 0785 WRITE (3,501) R 0786 501 FORMAT (1X,'VALUE OF SR/SV = ',G11.4,' USED TO CALCULATE CLDR IS 0787 1 OUTSIDE PREFERRED RANGE OF 0.0 TO 0.7') 0788 502 T=21.7949*R**5-46.4744*R**4+36.9347*R**3-14.259*R*R+3.70551*R-.000 0789 1057815 0790 CLDR=CLAVT*T*SV*ZV/(S*B) 0791 WRITE(*,7014) CLDA 7014 FORMAT(1X,'CLDA=',F15.8) GO TO 914 0792 C 0793 1015 R=SR/SV 0794 IF (R.GE.0.0.AND.R.LE..7) GO TO 512 0795 WRITE (3,511) R 0796 511 FORMAT (1X,'VALUE OF SR/SV = ',G11.4,' USED TO CALCULATE CNDR IS 0797 1 OUTSIDE PREFERRED RANGE OF 0.0 TO 0.7') 0798 512 T=21.7949*R**5-46.4744*R**4+36.9347*R**3-14.259*R*R+3.70551*R-.000 0799 1057815 0800 CNDR=-CLAVT*T*SV*LT*ETAV/(S*B) 0801 WRITE(*,7015) CNDR 7015 FORMAT(1X,'CNDR=',F15.8) GO TO 915 0802 C 0803 END 0804 C*********************************************************************** 0805 C 0806 C----------------------------------------------------------------------- C SUBROUTINE LATPRG 0807 C----------------------------------------------------------------------- 0808 C 0809 C THIS SUBROUTINE PREPARES THE DIMENSIONAL LATERAL 0810 C STABILITY DERIVATIVES FOR SUBMISSION TO THE SOLVE 0811 C ROUTINE OR THE ROOT LOCUS ROUTINE. 0812 C IF THESE ROUTINES ARE NOT AVAILABLE, DUMMIES 0813 C MAY BE SUBSTITUTED AND THE TIME RESPONSE AND BODE 0814 C PLOTS OBTAINED FROM LATOLD. 0815 C 0816 C DATE OF CURRENT VERSION: AUGUST 1981 0817 C 0818 C 0819 C*********************************************************************** 0820 C 0821 SUBROUTINE LATPRG(YV,U,GCOSGM,GSINGM,YP,YR,LB,LP,LR, 0822 1IXX,IZZ,IXZ,NB,NP,NR,YIN,LIN,NIN,IF, 0823 2AMP,ISET,W,W0,LK,LBB) 0824 IMPLICIT REAL*8(A-H,O-Z) 0825 REAL NG,DG,GAIND,W0,W,AMP 0826 REAL*8 NB,NP,NR,LB,LP,LR,IXX,IZZ,IXZ,LIN,NIN 0827 COMMON/TRANS/ JPT,JEND,IVERY 0828 DIMENSION A(3,3,4),B(3,4),DG(10),GN(3,10),NG(10),MZ(3) 0829 C 0830 ICNG=0 0831 ISW=0 0832 JWRITE=3 0833 C 0834 C ZERO OUT A-ARRAY 0835 C 0836 DO 1 I=1,3 0837 DO 1 J=1,3 0838 DO 1 K=1,4 0839 1 A(I,J,K)=0.0D0 0840 C 0841 C LOAD A-ARRAY 0842 C 0843 A(1,1,1)=-YV 0844 A(1,1,2)=1.0D0 0845 A(1,2,1)=-GCOSGM/U 0846 A(1,2,2)=-YP/U 0847 A(1,3,1)=GSINGM/U 0848 A(1,3,2)=1.0D0-YR/U 0849 A(2,1,1)=-LB 0850 A(2,2,2)=-LP 0851 A(2,2,3)=1.0D0 0852 A(2,3,2)=-LR 0853 A(2,3,3)=-IXZ/IXX 0854 A(3,1,1)=-NB 0855 A(3,2,2)=-NP 0856 A(3,2,3)=-IXZ/IZZ 0857 A(3,3,2)=-NR 0858 A(3,3,3)=1.0D0 0859 C 0860 C ZERO OUT B-VECTOR 0861 C 0862 DO 2 I=1,3 0863 DO 2 K=1,4 0864 2 B(I,K)=0.0D0 0865 C 0866 C LOAD B-VECTOR 0867 C 0868 B(1,1)=YIN/U 0869 B(2,1)=LIN 0870 B(3,1)=NIN 0871 IF(LBB.EQ.3) GO TO 6 0872 WRITE(JWRITE,3) 0873 WRITE(JWRITE,4) 0874 WRITE(JWRITE,5) 0875 3 FORMAT(///,10X,'THE FIRST VARIABLE IN THE SUBSEQUENT ', 0876 1'COMPUTATIONS IS THE SIDESLIP ANGLE') 0877 4 FORMAT(/,10X,'THE SECOND VARIABLE REPRESENTS THE BANK', 0878 1' ANGLE') 0879 5 FORMAT(/,10X,'THE THIRD VARIABLE REPRESENTS THE YAW ' 0880 1'ANGLE',///) 0881 6 CONTINUE 0882 GAIND=1.0 0883 C 0884 IF(LBB.GT.2) WRITE(JWRITE,7) LBB 0885 7 FORMAT(124X,'LBB=',I2) 0886 C 0887 C "SOLVE" IS DOCUMENTED IN THE BOOK "FORTRAN CODES 0888 C FOR CLASSICAL METHODS IN LINEAR DYNAMICS". 0889 C 0890 IF(LBB.LT.3) CALL SOLVE(A,B,GN,DG,ICNG,ISW,LK,MZ,MM,LBB, 0891 1ISET,IF,AMP,GAIND,W,W0) 0892 C 0893 IF(JPT.EQ.2) WRITE(JWRITE,8) ISW,ICNG,LK 0894 8 FORMAT(108X,'ISW=',I1,2X,'ICNG=',I1,2X,'LK=',I1) 0895 IF(LBB.GE.3) CALL CHANGE(A,B,GN,DG,ICNG,ISW,LK,MZ,MM,LBB, 0896 2ISET,IF,AMP,GAIND,W,W0) 0897 C 0898 RETURN 0899 END 0900 C************************************************************************ 0928 C C 0929 C______________________________________________________________________ C 0930 C SUBROUTINE CHANGE C 0931 C-----------------------------------------------------------------------C 0932 C C 0933 C THIS SUBROUTINE PREPARES THE DIMENSIONAL LONGITUDINAL C 0934 C STABILITY DERIVATIVES FOR USE IN THE LOCUS ROUTINE. C 0935 C C 0936 C C 0937 C NOTE: CRAMER, LOCUS, AND PPLOT ARE DOCUMENTED IN C 0938 C THE BOOK "FORTRAN CODES FOR CLASSICAL METHODS IN C 0939 C LINEAR DYNAMICS". C 0940 C C 0941 C C 0942 C DATE OF CURRENT VERSION: AUGUST 1981 C 0943 C C 0944 C C 0945 C***********************************************************************C 0946 C C 0947 SUBROUTINE CHANGE(A,B,GN,DG,ICNG,ISW,LK,MZ,MM,LBB,ISET,IF,AMP,GAIN 0948 1D,W,W0) 0949 IMPLICIT REAL*8(A-H,O-Z) 0950 REAL NG,DG,W0,W,AMP 0951 COMMON/TRANS/ JPT,JEND,IVERY 0952 COMMON/LOCO/ RR(200,9),RI(200,9) 0953 DIMENSION AA(10),AK(10),Z(20),ROOTR(10),ROOTI(10) 0954 DIMENSION A(3,3,4),B(3,4),X(10), DG(10),GN(3,10), 0955 1NG(10),MZ(3),UC(10),RRR(23,9),RRI(23,9) 0956 C 0957 JWRITE=3 0958 C 0959 IF(JPT.LT.4) WRITE(JWRITE,26) JPT,JEND,IVERY,LBB 0960 26 FORMAT(/,90X,'JPT=',I2,2X,'JEND=',I2,2X,'IVERY=',I2,2X, 0961 1'LBB=',I2,/) 0962 C 0963 CALL CRAMER(A,B,GN,DG,ICNG,ISW,LK,MZ,MM) 0964 C 0965 IF(LBB.EQ.5) RETURN 0966 C 0967 DO 27 I=1,MM 0968 AK(I)=0.0D0 0969 27 AA(I)=DG(I) 0970 DK=0.0D0 0971 ISW=4 0972 C 0973 CALL LOCUS(AA,AK,Z,DK,1,MM,MM1,MM2,JK,ISW) 0974 C 0975 DO 30 I=1,MM1 0976 MG=2*I-1 0977 ROOTR(I)=Z(MG) 0978 ROOTI(I)=Z(MG+1) 0979 30 CONTINUE 0980 DO 31 I=1,MM1 0981 RRR(JPT,I)=ROOTR(I) 0982 31 RRI(JPT,I)=ROOTI(I) 0983 C 0984 IF(JEND.EQ.0) RETURN 0985 C 0986 C LOAD COMMON BLOCK FOR PLOT ROUTINE. 0987 C 0988 DO 33 J=1,JPT 0989 DO 33 JJ=1,MM1 0990 RR(J,JJ)=RRR(J,JJ) 0991 33 RI(J,JJ)=RRI(J,JJ) 0992 WRITE(3,32) 0993 32 FORMAT(1H1) 0994 C 0995 CALL PPLOT(JPT,MM1) 0996 C 0997 IF(JPT.GT.22) CALL EXIT 0998 RETURN 0999 END 1000 C***********************************************************************C 0001 C C 0002 C-----------------------------------------------------------------------C 0003 C SUBROUTINE FORCEL C 0004 C-----------------------------------------------------------------------C 0005 C C 0006 C PURPOSE: C 0007 C C 0008 C TO DETERMINE AN AIRCRAFT'S STATIC LATERAL DIRECTIONAL C 0009 C STABILITY AND CONTROL FORCE CHARACTERISTICS GIVEN C 0010 C THE AIRCRAFT'S GEOMETRIC AND INERTIAL CHARACTERISTICS C 0011 C C 0012 C C 0013 C VARIABLE IDENTIFICATION: C 0014 C C 0015 C CAL - AILERON CHORD C 0016 C FEI - BANK ANGLE C 0017 C FA - AILERON STICK FORCE C 0018 C TI - TIME TO REACH A SPECIFIED BANK ANGLE C 0019 C YI - INBOARD EXTENT OF AILERON MEASURED ALONG WING C 0020 C SPAN FROM AIRCRAFT CENTERLINE. C 0021 C YP - SPAN-WISE LOCATION OF AILERON CENTER C 0022 C DELTR - RUDDER TAB DEFLECTION(DEGREES) C 0023 C G - AILERON GEAR RATIO C 0024 C GR - RUDDER GEAR RATIO C 0025 C BR - RUDDER SPAN C 0026 C FR - PEDAL FORCE C 0027 C FRDB - CHANGE IN PEDAL FORCE WITH CHANGE IN SIDESLIP C 0028 C ANGLE. C 0029 C DBDR - CHANGE IN SIDESLIP ANGLE WITH CHANGE IN C 0030 C RUDDER DEFLECTION C 0031 C BETA - SIDESLIP ANGLE(DEGREES) C 0032 C DELR - RUDDER DEFLECTION(DEGREES) C 0033 C DELA - AILERON DEFLECTION(DEGREES) C 0034 C IRUD = 1,CALCULATE AILERON-RELATED ITEMS ONLY C 0035 C = 2,CALCULATE RUDDER-RELATED ITEMS ONLY C 0036 C = 3,CALCULATE BOTH RUDDER AND AILERON ITEMS C 0037 C 0038 C-----------------------------------------------------------------------C 0039 C C 0040 C F1,F2,CA,CB,CB1,T,PHI,CCT,BBT,IGAP,CHALP,CHDEL, C 0041 C FC1,FC2,BE,BB,AR,IBAL,CHDELT,IREAD,IREAD1,DELCHA, C 0042 C AND DELCHD C 0043 C - HAVE THE SAME MEANING AS IN SUBROUTINE FORCE C 0044 C EXCEPT THAT THE CONTROL SURFACE DESCRIBED C 0045 C IS AN AILERON OR A RUDDER. C 0046 C C 0047 C-----------------------------------------------------------------------C 0048 C 0049 C CNB - YAWING MOMENT VARIATION WITH SIDESLIP ANGLE C 0050 C CNDA - YAWING MOMENT VARIATION WITH AILERON DEFLECTION C 0051 C CNDR - YAWING MOMENT VARIATION WITH RUDDER DEFLECTION C 0052 C CLB - ROLLING MOMENT VARIATION WITH SIDESLIP ANGLE C 0053 C CLP - ROLLING MOMENT VARIATION WITH ROLL RATE C 0054 C CLDA - ROLLING MOMENT VARIATION WITH AILERON DEFLECTION C 0055 C CLDR - ROLLING MOMENT VARIATION WITH RUDDER DEFLECTION C 0056 C C 0057 C OTHER PARAMETERS HAVE THE SAME DEFINITIONS AS IN C 0058 C SUBROUTINE LATDER C 0059 C C 0060 C SEE SUBROUTINE FORCE FOR OTHER PERTINANT COMMENTS. C 0061 C C 0062 C NOTES: C 0063 C C 0064 C 1. RIGHT PEDAL FORCE IS POSITIVE C 0065 C 2. NOSE LEFT SIDESLIP IS POSITIVE C 0066 C 3. RUDDER TRAILING EDGE LEFT IS POSITIVE C 0067 C 4. RIGHT WING DOWN IS A POSITIVE BANK ANGLE C 0068 C 5. LEFT AILERON TRAILING EDGE DOWN IS POSITIVE AILERON C 0069 C DEFLECTION C 0070 C 6. STICK RIGHT IS POSITIVE STICK FORCE. C 0071 C 7. NOSE RIGHT YAW IS POSITIVE C 0072 C 8. POSITIVE YAWING MOMENT IS NOSE RIGHT C 0073 C 9. POSITIVE ROLLING MOMENT IS RIGHT WING DOWN C 0074 C C 0075 C C 0076 C ALL COMPUTATIONS ARE CARRIED OUT IN DOUBLE PRECISION. C 0077 C C 0078 C DATE OF CURRENT VERSION: FEBRUARY 1982 C 0079 C C 0080 C C 0081 C THIS VERSION HAS BEEN MODIFIED TO COMPILE UNDER C MICROSOFT FORTRAN77 VERSION 3.20. CHANGES NECESSARY C TO PERMIT THIS ARE GENERALLY INDICATED BY THE C ABSENCE OF LINE NUMBERS. C C C DATE OF MODIFICATION: JUNE 1984 C C C***********************************************************************C 0082 C SUBROUTINE FORCEL(ETAV,SV,BV,U,RHO,MS,B,CH,BA,SR,CAL,YI, 0084 1CNB,CNDA,CLDR,CLDA,CNDR,CLB,IXX,CLP,S,LBB) 0085 IMPLICIT REAL*8(A-H,O-Z) 0086 REAL*8 MS,IXX 0087 C 0088 C1----------------------------------------------------------------------C 0089 READ(1,1) DELA,DELR,TI,G,DELTR,GR,BR 0090 1 FORMAT(7F10.3) 0091 C2----------------------------------------------------------------------C 0092 READ(1,2) IRUD,IREAD,IREAD1,BETA,CHDELT 0093 2 FORMAT(3I10,2F10.3) 0094 C-----------------------------------------------------------------------C 0095 C 0096 WRITE(3,10) 0097 10 FORMAT(1H1,////////,30X,'STATIC LATERAL-DIRECTIONAL ', 0098 1'STABILITY AND CONTROL CHARACTERISTICS',/,30X,64('-'),////) 0099 GO TO(21,22,21),IRUD 0100 C 0101 C3----------------------------------------------------------------------C 0102 21 READ(1,3) CA,CB,CB1,T,PHI,IREAD,IREAD1 0103 3 FORMAT(5F10.3,2I10) 0104 C4----------------------------------------------------------------------C 0105 READ(1,27) IBAL,IGAP 0106 27 FORMAT(2I10) 0107 C-----------------------------------------------------------------------C 0108 FACTOR=CAL/(CA+CB) 0109 IF((CA+CB).EQ.CAL) GO TO 11 0110 CA=CA*FACTOR 0111 CB=CB*FACTOR 0112 CB1=CB1*FACTOR 0113 T=T*FACTOR 0114 11 CONTINUE 0115 C5----------------------------------------------------------------------C 0116 IF(IREAD.EQ.1) READ(1,8) CHALP,CHDEL 0117 8 FORMAT(2F10.7) 0118 C-----------------------------------------------------------------------C 0119 IF(IREAD.EQ.1) GO TO 30 0120 FC1=0.4212D0 0121 FC2=-0.016848D0 0122 C6----------------------------------------------------------------------C 0123 IF(IREAD1.EQ.1) READ(1,4) FC1,FC2 0124 4 FORMAT(2F10.3) 0125 C-----------------------------------------------------------------------C 0126 F2=1.0D0-DSQRT(1.0D0-((1.0D0+CB1/CA)/(1.0D0+CB/CA))**2) 0127 F1F2=(FC1*(CB/CA)**2+FC2)*(0.4D0*CA/T)**0.58D0 0128 F1=F1F2/F2 0129 C7----------------------------------------------------------------------C 0130 IF(IBAL.NE.1) READ(1,4) BE,BB 0131 C-----------------------------------------------------------------------C 0132 IF(IBAL.NE.1) F1=F1*BB/BE 0133 IF(IGAP.EQ.0) GO TO 5 0134 DELCHA=0.017D0*F1+5.0D-6*PHI 0135 DELCHD=0.10D0*F1F2+1.7D-4*PHI 0136 GO TO 6 0137 5 DELCHA=0.14D0*(CA/CH)**2*F1 0138 DELCHD=0.09D0*DSQRT(DABS(CA/CH))*F1 0139 6 CONTINUE 0140 C-----------------------------------------------------------------------C 0141 C 0142 C CALCULATE HINGE MOMENT COEFFICIENTS 0143 C 0144 C-----------------------------------------------------------------------C 0145 AR=B*B/S 0146 CHALP=AR/(AR+2.0D0)*(5.0D-4*PHI-0.02D0*(CA/CH)+DELCHA 0147 1-0.006D0) 0148 CHDEL=AR/(AR+2.0D0)*(4.0D-4*PHI-0.02D0*(CA/CH)+DELCHD 0149 1-0.0106D0) 0150 30 CONTINUE 0151 IF(CHDEL.GT.0.0D0) WRITE(3,12) CHDEL,DELCHD 0152 12 FORMAT(//,20X,'PROBABLE USER ERROR: CHDEL=',F17.7,3X, 0153 1'GREATER THAN ZERO',3X,'DELCHD=',F17.7,//) 0154 C-----------------------------------------------------------------------C 0155 C 0156 C CALCULATE BANK ANGLE ACHIEVED AND AILERON FORCE REQUIRED 0157 C 0158 C-----------------------------------------------------------------------C 0159 ARGU=CLP*RHO*U*S*S/4.0D0/IXX 0160 IF(ARGU*TI.GT.63.0D0) ARGU=63.0D0/TI 0161 IF(ARGU*TI.LT.-63.0D0) ARGU=-63.0D0/TI 0162 IF(DABS(ARGU*TI).LT.1.0D-63) ARGU=0.0D0 0163 FEI=-CLDA/CLP*U/B*DELA*TI-(1.0D0-DEXP(ARGU*TI))/ARGU/CLP*U/B 0164 1*CLDA*2.0D0*DELA 0165 YP=YI+BA/2.0D0 0166 EN=FEI/TI*YP/U/DELA 0167 FA=-RHO*BA*CAL*CAL*U*U/2.0D0*G*DELA*(1.0D0-2.0D0*EN*CHALP/ 0168 1CHDEL)*CHDEL 0169 C-----------------------------------------------------------------------C 0170 C 0171 C WRITE RESULTS FOR AILERON INPUTS 0172 C 0173 C-----------------------------------------------------------------------C 0174 WRITE(3,13) FEI,TI,FA 0175 13 FORMAT(10X,112('*'),/,10X,'*',110(' '),'*',/,10X,'*',10X,'BANK ', 0176 1'ANGLE ACHIEVED =',F17.7,3X,'DEGREES',3X,'AT TIME =',F17.7,3X,'SEC 0177 1ONDS',13X,'*',/,10X,'*',10X,'AILERON STICK FORCE REQUIRED DURING 0178 1ROLL =',F17.7,40X,'*',/,10X,'*',110(' '),'*',/,10X,112('*')) 0179 WRITE(3,28) 0180 WRITE(3,29) CHALP,CHDEL,DELA,DELR,TI,G,DELTR,GR,BR 0181 28 FORMAT(////,35X,'COMPUTED AND READ IN CHARACTERISTICS',/,35X, 0182 136('-')) 0183 29 FORMAT(/,35X,'CHALP =',27X,F17.7,/,35X,'CHDEL =',27X,F17.7, 0184 2,/,35X,'AILERON DEFLECTION =',14X,F17.7,/,35X,'RUDDER', 0185 3' DEFLECTION =',15X,F17.7,/,35X,'TIME TO REACH A GIVEN BANK', 0186 4' ANGLE =',F17.7,/,35X,'AILERON GEAR RATIO =',14X,F17.7,/, 0187 535X,'RUDDER TAB DEFLECTION =',11X,F17.7,/,35X,'RUDDER', 0188 6' GEAR RATIO =',15X,F17.7,/,35X,'RUDDER SPAN =',21X,F17.7) 0189 WRITE(3,15) CA,CB,CB1,T,PHI,IGAP,IBAL 0190 15 FORMAT(9X,'AILERON ATTIRBUTES: CA=',F5.2,3X,'CB=',F5.2,3X, 0191 1'CB1=',F5.2,3X,'T=',F5.2,3X,'PHI=',F5.2,3X,'IGAP=',I2, 0192 23X,'IBAL=',I3) 0193 IF(F1F2.LT.0.0D0) WRITE(3,14) 0194 14 FORMAT(/,25X,'WARNING: SPECIFIED GEOMETRIC PARAMETERS', 0195 1' YIELD NEGATIVE VALUE FOR F1F2. INCREASE CB.',/) 0196 IF(IRUD.LT.3) GO TO 24 0197 C 0198 C-----------------------------------------------------------------------C 0199 C 0200 C BEGIN COMPUTATIONS FOR RUDDER FORCES AND GRADIENTS 0201 C 0202 C8----------------------------------------------------------------------C 0203 22 READ(1,3) CA,CB,CB1,T,PHI,IREAD,IREAD1 0204 C9----------------------------------------------------------------------C 0205 READ(1,27) IBAL,IGAP 0206 C-----------------------------------------------------------------------C 0207 CR=SR/BR 0208 FACTOR=CR/(CA+CB) 0209 IF((CA+CB).EQ.CR) GO TO 41 0210 CA=CA*FACTOR 0211 CB=CB*FACTOR 0212 CB1=CB1*FACTOR 0213 T=T*FACTOR 0214 41 CONTINUE 0215 CV=SV/BV 0216 AR=BV*BV/SV 0217 C10---------------------------------------------------------------------C 0218 IF(IREAD.EQ.1) READ(1,8) CHALP,CHDEL 0219 C-----------------------------------------------------------------------C 0220 IF(IREAD.EQ.1) GO TO 33 0221 FC1=0.4212D0 0222 FC2=-0.016848D0 0223 C11---------------------------------------------------------------------C 0224 IF(IREAD1.EQ.1) READ(1,4) FC1,FC2 0225 C-----------------------------------------------------------------------C 0226 C 0227 F2=1.0D0-DSQRT(1.0D0-((1.0D0+CB1/CA)/(1.0D0+CB/CA))**2) 0228 F1F2=(FC1*(CB/CA)**2+FC2)*(0.4D0*CA/T)**0.58D0 0229 F1=F1F2/F2 0230 C 0231 C12---------------------------------------------------------------------C 0232 IF(IBAL.NE.1) READ(1,4) BE,BB 0233 C-----------------------------------------------------------------------C 0234 C 0235 IF(IBAL.NE.1) F1=F1*BB/BE 0236 IF(IGAP.EQ.0) GO TO 51 0237 DELCHA=0.017D0*F1+5.0D-6*PHI 0238 DELCHD=0.10D0*F1F2+1.7D-4*PHI 0239 GO TO 61 0240 51 DELCHA=0.14D0*(CA/CV)**2*F1 0241 DELCHD=0.09D0*DSQRT(DABS(CA/CV))*F1 0242 61 CONTINUE 0243 C-----------------------------------------------------------------------C 0244 C 0245 C CALCULATE HINGE MOMENT COEFFICIENTS 0246 C 0247 C-----------------------------------------------------------------------C 0248 CHALP=AR/(AR+2.0D0)*(5.0D-4*PHI-0.02D0*(CA/CV)+DELCHA 0249 1-0.006D0) 0250 CHDEL=AR/(AR+2.0D0)*(4.0D-4*PHI-0.02D0*(CA/CV)+DELCHD 0251 1-0.0106D0) 0252 33 CONTINUE 0253 IF(CHDELT.NE.0.0D0) GO TO 23 0254 C13---------------------------------------------------------------------C 0255 READ(1,4) CCT,BBT 0256 C-----------------------------------------------------------------------C 0257 C 0258 CHDELT=(0.004D0*(CCT/CR)-0.000833D0*(CR/CV)-0.014534D0) 0259 1*BBT/BR 0260 23 CONTINUE 0261 IF(CHDEL.GT.0.0D0) WRITE(3,12) CHDEL,DELCHD 0262 C-----------------------------------------------------------------------C 0263 C 0264 C DETERMINE RUDDER FORCE, CHANGE IN RUDDER FORCE WITH CHANGE 0265 C IN SIDESLIP, AND CHANGE IN SIDESLIP WITH RUDDER DEFLECTION 0266 C 0267 C-----------------------------------------------------------------------C 0268 CR=SR/BR 0269 FR=GR*RHO*U*U/2.0D0*ETAV*SR*CR*(-CHALP*BETA+CHDEL*DELR 0270 1+CHDELT*DELTR) 0271 FRDB=-GR*RHO*U*U/2.0D0*ETAV*SR*CR*(CHALP+CHDEL*CNB/CNDR) 0272 DBDR=-CNDR/CNB 0273 C-----------------------------------------------------------------------C 0274 C 0275 C WRITE RESULTS FOR RUDDER INPUTS 0276 C 0277 C-----------------------------------------------------------------------C 0278 WRITE(3,7) FR,BETA,FRDB,DBDR 0279 7 FORMAT(////,10X,112('*'),/,10X,'*',10X,'RUDDER FORCE', 0280 1' REQUIRED =',F17.7,3X,'AT SIDESLIP ANGLE =',F17.7,21X,'*',,/, 0281 210X,'*',10X,'CHANGE IN RUDDER FORCE WITH CHANGE IN SIDESLIP', 0282 3' ANGLE =',F17.7,29X,'*',/,10X,'*',10X,'CHANGE IN SIDESLIP', 0283 4' ANGLE WITH CHANGE IN RUDDER DEFLECTION =',F17.7,24X,'*', 0284 5/,10X,112('*')) 0285 WRITE(3,9)CHALP,CHDEL,CHDELT 0286 9 FORMAT(////,20X, 'COMPUTED RESULTS FOR RUDDER HINGE MOMENTS',/, 0287 120X,41('-'),/,20X,'CHALP =',F17.7,/,20X,'CHDEL =',F17.7,/, 0288 220X,'CHDELT =',F17.7) 0289 WRITE(3,16) CA,CB,CB1,T,PHI,IGAP,IBAL 0290 16 FORMAT(9X,'RUDDER ATTIRBUTES: CA=',F5.2,3X,'CB=',F5.2,3X, 0291 1'CB1=',F5.2,3X,'T=',F5.2,3X,'PHI=',F5.2,3X,'IGAP=',I2, 0292 23X,'IBAL=',I3) 0293 IF(F1F2.LT.0.0D0) WRITE(3,14) 0294 C 0295 24 IF(LBB.EQ.6) STOP 0296 IF(LBB.EQ.7) LBB=2 0297 IF(LBB.EQ.8) LBB=3 0298 RETURN 0299 END 0300 C***********************************************************************C C***********************************************************************C C************ ************C C******** ++++++++++++++++ ********C C****** ++++++++++++++++ ******C C**** ****C C** *** *** * * * **** **C C** * * * * * * * * **C C** * * * * * * * **C C** *** * * * * * **** **C C** * * * * * * * **C C** * * * * * * * * **C C** *** *** **** * **** **C C** **C C**** ****C C******* *******C C************ ************C C***********************************************************************C C***********************************************************************C C C C THIS VERSION IS A SUBSET OF SOLVE83 INTENDED C C SPECIFICALLY FOR USE WITH STADER AND LATDER. C C C C DATE OF PRESENT VERSION: JUNE 1984 C C C C***********************************************************************C C SUBROUTINE SOLVE(A,B,GN,DG,ICNG,ISW,LK,MZ,MM,LB,ISET,IF,AMP,GAIND, 1W,W0) IMPLICIT REAL*8(A-H,O-Z) REAL NG,DG,W0,W,AMP DIMENSION AA(10),AK(10),Z(20),ROOTR(10),ROOTI(10),NG(10) DIMENSION A(3,3,4),B(3,4),X(10),MZ(3),DG(10),GN(3,10),UC(10) C ISW=0 CALL CRAMER(A,B,GN,DG,ICNG,ISW,LK,MZ,MM) WRITE(*,1001) 1001 FORMAT(4X,'returned from CRAMER') IF(LB.EQ.0) GO TO 20 DO 18 JN=1,LK MZZ=MZ(JN) DO 17 L=1,MZZ 17 NG(L)=GN(JN,L) 18 CALL BODE(NG,DG,MZ(JN),MM) WRITE(*,1002) 1002 FORMAT(5X,'returned from BODE') 20 CONTINUE DO 27 I=1,MM AK(I)=0.0D0 27 AA(I)=DG(I) DK=0.0D0 WRITE(*,1003) 1003 FORMAT(6X,'calling subroutine LOCUS') CALL LOCUS(AA,AK,Z,DK,1,MM,MM1,MM2,JK,ISW) C DO 30 I=1,MM1 MG=2*I-1 ROOTR(I)=Z(MG) ROOTI(I)=Z(MG+1) 30 CONTINUE GAIND=DG(MM) DO 29 LG=1,LK MZZ=MZ(LG) DO 28 L=1,MZZ 28 UC(L)=GN(LG,L) 29 CALL RESVAX(ROOTR,ROOTI,AMP,GAIND,IF,UC,MM1,MZZ,W,ISET,W0) RETURN END C***********************************************************************C C C C PROGRAM CRAMER C C C C USES CRAMER'S RULE TO FIND TRANSFER FUNCTIONS FROM C C TRANSFORMED EQUATIONS OF MOTION C C***********************************************************************C C SUBROUTINE CRAMER(A,B,GN,DG,ICNG,ISW,LK,MZ,MM) IMPLICIT REAL*8(A-H,O-Z) REAL NG,DG DIMENSION A(3,3,4),B(3,4),X(10),Z(3,3,4),DG(10),GN(3,10), 1NG(10),MZ(3) C C EVALUATE DENOMINATOR DETERMINANT C WRITE(*,1000) 1000 FORMAT(3X,'subroutine CRAMER invoked') CALL DETERM(A,X,ISW) M=10 3 FORMAT(4F10.7) DO 5 N=1,10 M=11-N IF(X(M).NE.0.0D0) GO TO 6 5 CONTINUE 6 DO 7 I=1,M 7 DG(I)=X(I) DNORM=DG(M) DO 22 I=1,M 22 DG(I)=DG(I)/DNORM MM=M DO 13 L=1,LK DO 8 I=1,3 DO 8 J=1,3 DO 8 K=1,4 8 Z(I,J,K)=A(I,J,K) DO 9 I=1,3 DO 9 K=1,4 9 Z(I,L,K)=B(I,K) C C EVALUATE NUMERATOR DETERMINANTS C CALL DETERM(Z,X,ISW) M=10 DO 10 N=1,10 M=11-N IF(X(M).NE.0.0D0) GO TO 11 10 CONTINUE 11 DO 12 I=1,M 12 GN(L,I)=X(I)/DNORM MZ(L)=M 13 CONTINUE DO 15 L=1,LK MZZ=MZ(L) WRITE(3,21) WRITE(3,14) (GN(L,I),I=1,MZZ) WRITE(3,19) L WRITE(3,14) (DG(I),I=1,MM) 15 WRITE(3,21) 14 FORMAT( 10X,F10.5,'+',F10.5,'*S+',F10.5,'*S**2+',F10.5, 1'*S**3+',F10.5,'*S**4+',F10.5,'*S**5+',F10.5,'*S**6+', 2F10.5,'*S**7+',F10.5,'*S**8') 19 FORMAT(4X,'TF(',I1,')=',110('-')) 21 FORMAT(//) RETURN END C***********************************************************************C C C C PROGRAM DETERM C C C C***********************************************************************C C SUBROUTINE DETERM(E,X,ISW) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(3,4),Z(10),X(10),E(3,3,4),IY(18) DATA IY/1,2,3,2,3,1,3,1,2,3,2,1,1,3,2,2,1,3/ C M=10 DO 69 N=1,10 69 X(N)=0.0D0 DO 70 K=3,18,3 DO 73 I=1,3 DO 73 J=1,4 L=K-3+I 73 A(I,J)=E(I,IY(L),J) CALL MUL(A,Z,M,ISW) DO 71 N=1,10 IF(L.GE.12) GO TO 72 X(N)=X(N)+Z(N) GO TO 71 72 X(N)=X(N)-Z(N) 71 CONTINUE 70 CONTINUE DO 80 N=1,10 M=11-N IF(X(M).NE.0.0D0) GO TO 90 80 CONTINUE 90 CONTINUE RETURN END C***********************************************************************C C C C PROGRAM MUL C C C C***********************************************************************C C SUBROUTINE MUL(A,Z,M,ISW) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(3,4),Z(10) C Z(1)=A(1,4)*A(2,4) Z(2)=A(2,4)*A(1,3)+A(1,4)*A(2,3) Z(3)=A(1,2)*A(2,4)+A(1,3)*A(2,3)+A(2,2)*A(1,4) Z(4)=A(2,4)*A(1,1)+A(1,2)*A(2,3)+A(1,3)*A(2,2)+A(1,4)*A(2,1) Z(5)=A(1,1)*A(2,3)+A(1,2)*A(2,2)+A(1,3)*A(2,1) Z(6)=A(2,2)*A(1,1)+A(1,2)*A(2,1) Z(7)=A(1,1)*A(2,1) X1=Z(7)*A(3,1) X2=A(3,4)*Z(1) X3=0.0D0 X4=0.0D0 X5=0.0D0 X6=0.0D0 X7=0.0D0 X8=0.0D0 X9=0.0D0 X10=0.0D0 DO 1 I=1,2 1 X3=X3+A(3,I)*Z(5+I) DO 2 I=1,3 2 X4=X4+A(3,I)*Z(I+4) DO 3 I=1,4 3 X5=X5+A(3,I)*Z(3+I) DO 4 I=1,4 4 X6=X6+A(3,I)*Z(2+I) DO 5 I=1,4 5 X7=X7+A(3,I)*Z(1+I) DO 6 I=1,4 6 X8=X8+A(3,I)*Z(I) DO 7 I=1,3 7 X9=X9+A(3,I+1)*Z(I) DO 8 I=1,2 8 X10=X10+A(3,I+2)*Z(I) Z(1)=X1 Z(2)=X3 Z(3)=X4 Z(4)=X5 Z(5)=X6 Z(6)=X7 Z(7)=X8 Z(8)=X9 Z(9)=X10 Z(10)=X2 RETURN END C***********************************************************************C C C C PROGRAM BODE C C C C EVALUATES AND PLOTS FREQUENCY RESPONSE (BODE) AMPLITUDE C C RATIO AND PHASE ANGLE C C***********************************************************************C C SUBROUTINE BODE(NG,DG,ING,IDG) IMPLICIT REAL*4 (N),COMPLEX*8 (C) INTEGER*4 NCYCLE COMMON/GOO/IK,IRUB,JZ,KJ DIMENSION NG(10),DG(10),DB(200),W(200),PHASE(200),XAXIS(200) C NCYCLE=5 C WRITE(*,1000) 1000 FORMAT(3X,'subroutine BODE invoked') CALL FREQ(W,IW,XAXIS,NCYCLE) WRITE(3,98) 98 FORMAT('1'////,10X,'THE COEFFICIENTS OF THE NUMERATOR'/) DO 70 I=1,ING 70 WRITE(3,101)I,NG(I) 101 FORMAT(10X,'NG(',I1,')=',F17.7) WRITE(3,102) 102 FORMAT(////10X,'THE COEFFICIENTS OF THE DENOMINATOR'/) DO 80 I=1,IDG 80 WRITE(3,104)I,DG(I) 104 FORMAT(10X,'DG(',I1,')='F17.7) CALL BODEC(W,IW,DB,PHASE,NG,ING,DG,IDG) WRITE(3,811) 811 FORMAT(1H1) WRITE(3,611) 611 FORMAT(////8X,'AXIS',14X,'W',16X,'DB',13X,'PHASE', 18X,'AXIS',14X,'W',16X,'DB',13X,'PHASE',/) DO 91 I=2,51 J=I+50 91 WRITE(3,601) XAXIS(I),W(I),DB(I),PHASE(I),XAXIS(J),W(J),DB(J), 1PHASE(J) 601 FORMAT(1X,8F16.7) IRUB=3 C C PLOT AMPLITUDE RATIO C CALL GENPT(XAXIS,DB,IW,0) IRUB=4 C C PLOT PHASE ANGLE C CALL GENPT(XAXIS,PHASE,IW,0) 3 RETURN END C***********************************************************************C C C C SUBROUTINE BODEC CALCULATES THE AMPLITUDE RATIO C C AND PHASE ANGLE. C C C C***********************************************************************C C SUBROUTINE BODEC(W,IW,DB,PHASE,NG,ING,DG,IDG) IMPLICIT REAL*4 (N),COMPLEX*8 (C) COMPLEX*8 JW REAL*8 JWW(2),DCN(2),DCD(2) REAL*4 DCNN(2),DCDD(2) REAL*8 DGG(10),GN(10),A,B DIMENSION NG(10),DG(10),DB(200),W(200),PHASE(200) C DO 1 J=1,10 DGG(J)=DG(J) 1 GN(J)=NG(J) DO 90 J=1,IW JW=CMPLX(0.0D0,W(J)) JWW(1)=0.0 JWW(2)=W(J) CALL CQPVAL(DCN,JWW,GN,ING) CALL CQPVAL(DCD,JWW,DGG,IDG) DCNN(1)=DCN(1) DCNN(2)=DCN(2) DCDD(1)=DCD(1) DCDD(2)=DCD(2) CN=CMPLX(DCNN(1),DCNN(2)) CD=CMPLX(DCDD(1),DCDD(2)) CG=CN/CD CALL HANGE(CG,A,B) GMAG=SQRT(A*A+B*B) DB(J)=20.*ALOG10(GMAG) PHASE(J)=ATAN2(B,A)*57.2958 IF((IDG.GT.(ING+2)).AND.B.GE.0.0) PHASE(J)=PHASE(J)-360.0 IF(IDG.GT.(ING+2)) GO TO 90 IF((PHASE(1).LE.-90.0).AND.B.GE.0.0) PHASE(J)=PHASE(J)-360.0 IF(J.LE.2) GO TO 90 IF(((PHASE(J-1)-PHASE(J-2)).GT.0.0).AND.(PHASE(1).GE.0.0 1.AND.B.LE.0.0)) PHASE(J)=PHASE(J)+360.0 IF((ABS(PHASE(J)-PHASE(J-1)).GT.170.0).AND.PHASE(J-1).GT. 1180.0) PHASE(J)=PHASE(J)+360.0 IF((ABS(PHASE(J)-PHASE(J-1)).GT.170.0).AND.PHASE(J-1).LT. 1-100.0) PHASE(J)=PHASE(J)-360.0 90 CONTINUE RETURN END C***********************************************************************C C C C SUBROUTINE HANGE C C C C THIS SUBROUTINE HAS BEEN INCLUDED BECAUSE THE MICROSOFT C C COMPILER INTRINSIC FUNCTIONS DREAL AND DIMAG CONTAIN AN C C ERROR. C C C C CURRENT VERSION: NOVEMBER 1988 C C C C***********************************************************************C C C SUBROUTINE HANGE(OUT1,RREAL,UNREAL) REAL*8 RREAL,UNREAL REAL*4 HOLD COMPLEX*8 OUT1,COMP DIMENSION HOLD(2) EQUIVALENCE (HOLD(1),COMP) C COMP=OUT1 RREAL=DBLE(HOLD(1)) UNREAL=DBLE(HOLD(2)) RETURN END SUBROUTINE FREQ(W,IW,XAXIS,NCYCLE) IMPLICIT REAL*4(A-H,O-Z) INTEGER*4 NCYCLE,NEXP DIMENSION W(200),F(20),XAXIS(200) DATA F/1.,1.1,1.25,1.4,1.6,1.8,2.,2.2,2.5,2.8,3.2,3.6,4.0,4.4,5., $ 5.6,6.4,7.,8.,9./ C DO 2 I=1,NCYCLE NEXP=I-4 WRITE(*,4) NEXP 4 FORMAT(1X,'NEXP=',I4) DO 1 J=1,20 K=20*I+J-20 W(K)=F(J)*10.**NEXP 1 XAXIS(K)=ALOG10(W(K)) 2 CONTINUE IW=NCYCLE*20+1 WRITE(*,5) IW 5 FORMAT(10X,'IW=',I4) W(IW)=(10.0)**(NEXP+1) XAXIS(IW)=ALOG10(W(IW)) RETURN END C********************************************************************** C C SUBROUTINE RESVAX C C PROVIDES COMPUTATION OF INVERSE LAPLACE TRANSFORM OF C TRANSFER FUNCTIONS WITHOUT USING COMPLEX ARTHMETIC C********************************************************************** SUBROUTINE RESVAX(ROOTR,ROOTI,AMP,GAIND,IF,NG,IDGM1,ING,W, 1ISET,W0) IMPLICIT REAL*8(N) REAL*8 ROOTR,ROOTI,GAIND DIMENSION NG(10),DG(10),ROOTR(10),ROOTI(10),OUTPUT(200),T(200),FOR 1CE(6,2) C WRITE(*,1999) 1999 FORMAT(1X,'program has entered subroutine RESVAX') SG=1.0 P2=6.28 W1=W0*P2 SR2=(10.0)**0.25 ICNT=1 IJKL=IDGM1 EPP=-0.000330 DO 7 I=1,IDGM1 7 IF(ROOTR(I).EQ.0.0.AND.ROOTI(I).EQ.0.0) ROOTR(I)=(7.*I-I*I-8.)*EPP 1/2.00159 DO 9 I=1,IDGM1 RIMAG=0.0 ISAG=0 IBAG=0 DO 4 J=ICNT,IDGM1 IF(I.EQ.J) GO TO 4 IF(ROOTR(J).EQ.0.0) GO TO 4 RREAL=ROOTR(I)/ROOTR(J) IF(ROOTI(J).EQ.0.0) RIMAG=0.001 RRIMAG=ROOTI(I)/(ROOTI(J)+RIMAG) IF(ABS(RREAL-1.0).LT..001) ISAG=1 IF(ABS(RRIMAG+1.).LT..001) IBAG=1 IF(ISAG.EQ.1.AND.IBAG.NE.1) ROOTR(I)=ROOTR(I)+EPP 4 CONTINUE 9 ICNT=ICNT+1 5 IF(ISET.EQ.1) GO TO 6 IF(IF.EQ.5) W=W1*SG IF(IF.EQ.4) W=1.0 IF(IF.EQ.3) W=2.0 6 CONTINUE CALL UNSPON(ROOTR,ROOTI,AMP,GAIND,IF,NG,IDGM1,ING,W) IF(IF.NE.5) GO TO 8 IDGM1=IDGM1-2 IF(W.GT.80.0) GO TO 8 IF(ISET.EQ.1) GO TO 8 SG=SG*SR2 GO TO 5 8 CONTINUE IDGM1=IJKL RETURN END C***********************************************************************C C C C SUBROUTINE UNSPON C C C C DETERMINES NUMBER OF PLOT SEGMENTS; PRINTS RESULTS C C***********************************************************************C C SUBROUTINE UNSPON(ROOTR,ROOTI,AMP,GAIND,IF,NG,IDGM1,ING,W) IMPLICIT REAL*8(N) REAL*8 ROOTR,ROOTI,GAIND INTEGER HILIM,NUM LOGICAL INRNG CHARACTER*4 FORCE DIMENSION NG(10),DG(10),ROOTR(10),ROOTI(10),OUTPUT(200),T(200),FOR 1CE(7,2) COMMON TDEL,TMAX,SMALL,IRE COMMON/GOO/IK,IRUB,JZ,KJ DATA FORCE/'IMPL','STEP','RAMP','PULS','RAMP','SINU','ACCE', 1'LSE',' ',' ' 1,'E','STEP','SOID','LERA'/ INRNG(LOWLIM,NUM,HILIM)=(NUM.GE.LOWLIM).AND.(NUM.LE.HILIM) C JZ=0 RTIME=2. PTIME=2. IERR=0 IF(W.EQ.0.0) GO TO 1 IF(IF.EQ.3) PTIME=W IF(IF.EQ.4) RTIME=W GO TO 4 1 IF(IF.EQ.1) IERR=1 IF(IF.EQ.2) IERR=2 IF(IF.EQ.6) IERR=3 4 IRE=0 SMALL=0.0 C C IF NO NUMERATOR COEFFICIENTS ARE PRESENT, TERMINATE C COMPUTATION. C IF(ING.LE.0) WRITE(3,64) 64 FORMAT(/,15X,'PROBABLE USER ERROR: NUMERATOR OF TRANSFER' 1' FUNCTION DOES NOT EXIST',/,30X,'COMPUTATION TERMINATING') IF(ING.LE.0) RETURN 5 IF(AMP.EQ.0.0) AMP=1.0 AM1=AMP/GAIND IZIP=0 IBANG=0 30 CONTINUE C CALL TIME(IF,AM1,PTIME,RTIME,W,OUTPUT,T,IO,NG,ING,ROOTR,ROOTI,IDGM 11,ICODE) C IF(ICODE.NE.0) GO TO 101 6 FORMAT ('1'////,10X,'THE COEFFICIENTS OF THE NUMERATOR'/) C C DETERMINE WHETHER ADDITIONAL PLOTS ARE NEEDED. C IF(TMAX*SMALL.LT.5.9) IRE=1 IF(TMAX*SMALL.GT.6.8) IRE=0 IF(IRE.GE.1) IZIP=IZIP+1 IF(IRE.EQ.0.AND.IZIP.GE.1) GO TO 23 IF(IRE.EQ.0) GO TO 31 IF(IZIP.GT.1) IRE=IZIP 31 CONTINUE IF(IRE.GE.2) GO TO 26 C C WRITE INPUT DATA C WRITE (3,6) DO 7 I=1,ING 7 WRITE (3,8) I,NG(I) 8 FORMAT (10X,'NG(',I1,')=',F17.7) WRITE (3,9) 9 FORMAT (////10X,'THE ROOTS OF THE DENOMINATOR'/) DO 10 I=1,IDGM1 10 WRITE (3,11) I,ROOTR(I),ROOTI(I) 11 FORMAT (10X,'ROOT(',I1,')=',F17.7,' + J ',F17.7) WRITE (3,12) IF,(FORCE(IF+1,I),I=1,2),AMP 12 FORMAT (////10X,'THE FORCING FUNCTION INDICATOR (IF) =',I5,//10X,' 1THIS IMPLIES THAT A ',2A4,' INPUT WAS USED.',//20X,'AMPLITUDE=',F1 17.7) IF(IF.NE.5) GO TO 26 WRITE(3,15) W WRITE(3,49) 49 FORMAT(1H1) 26 IF(IF.EQ.5.AND.IRE.GE.2) WRITE(3,33) W IF(IRE.GE.2) GO TO 51 33 FORMAT(1H1,10X,'THE SYSTEM IS BEING FORCED BY A SINE WAVE', 1' OF FREQUENCY=',F17.7) IF (IF.EQ.3) WRITE (3,13) PTIME 13 FORMAT (/20X,'PTIME=',F17.7) IF (IF.EQ.4) WRITE (3,14) RTIME IF(INRNG(1,IF,2).AND.INRNG(1,IERR,2)) WRITE(3,52) IF(IF.EQ.6.AND.IERR.EQ.3) WRITE(3,52) IF(INRNG(3,IF,4)) WRITE(3,49) 52 FORMAT(//,20X,'THE OUTPUT HAS BEEN SPECIFIED TO BE THE', 1' DIFFERENCE BETWEEN THE RESPONSE AND THE COMMAND, IE., THE ERROR' 2,//) 14 FORMAT (/20X,'RTIME WAS ADJUSTED TO ',F10.4,' SECONDS.') 15 FORMAT (/20X,'FREQUENCY=',F17.7) 51 CONTINUE C C ADJUST ROOT COUNTERS IF PREMATURE TERMINATION OCCURS C IN TYME. C IF(IRE.GE.2.AND.INRNG(3,IF,4)) WRITE(3,49) 101 IF(ICODE.NE.2) GO TO 100 WRITE(3,103) IDGM1 WRITE(3,102) (ROOTR(LL),LL=1,IDGM1) WRITE(3,102) (ROOTI(LL),LL=1,IDGM1) 102 FORMAT(3X,5(E15.8,2X)) 103 FORMAT(50X,I2) IF(IF.EQ.0) GO TO 77 GO TO(71,72,73,74,75,76),IF 71 IDGM1=IDGM1-2 GO TO 77 72 IDGM1=IDGM1-3 GO TO 77 73 IDGM1=IDGM1-2 GO TO 77 74 IDGM1=IDGM1-3 GO TO 77 75 IDGM1=IDGM1-3 GO TO 77 76 IDGM1=IDGM1-4 77 CONTINUE WRITE(3,102) (ROOTR(LL),LL=1,IDGM1) WRITE(3,102) (ROOTI(LL),LL=1,IDGM1) IBANG=IBANG+1 IF(IBANG.GT.3) GO TO 100 GO TO 30 100 CONTINUE IF(ICODE.EQ.0) GO TO 83 WRITE(3,16) ICODE 16 FORMAT (//,' ICODE =',I10) IF(ICODE-2) 80,81,82 80 WRITE(3,17) 17 FORMAT (//10X,'THE COMPLEX PART OF THE OUTPUT VECTOR BECAME 1 SIGNIFICANT') RETURN 81 WRITE(3,18) 18 FORMAT (//10X,'MULTIPLE ROOTS ENCOUNTERED') RETURN 82 WRITE(3,19) 19 FORMAT (//10X,'BAD ENTRY - CHECK POLYNOMIAL ORDERS OF T.F.') RETURN 83 IF(IF.LT.3) WRITE(3,49) C C COMPUTE SYSTEM RESPONSE ERROR, IF IERR CODE IS SET. C IF(IERR.EQ.0) GO TO 44 IF(IERR.EQ.3) GO TO 45 IF(IERR.EQ.2) GO TO 42 C C CHANGE OUTPUT TO ERROR RESPONSE IF CODES ARE SET. C DO 43 I=2,IO 43 OUTPUT(I)=OUTPUT(I)-AM1 GO TO 44 45 DO 46 I=1,IO 46 OUTPUT(I)=OUTPUT(I)-0.5*AM1*T(I)*T(I) GO TO 44 42 DO 41 I=1,IO 41 OUTPUT(I)=OUTPUT(I)-AM1*T(I) 44 CONTINUE C C WRITE OUTPUT C WRITE (3,20) JF=0 JA=IO/4.0 JB=2*JA JC=3*JA DO 21 I=1,JA JF=JC+I 21 WRITE(3,22) T(I),OUTPUT(I),T(JA+I),OUTPUT(JA+I),T(JB+I), 1OUTPUT(JB+I),T(JC+I),OUTPUT(JC+I) 20 FORMAT(//,8X,'TIME',13X,'OUTPUT',8X,'TIME',13X,'OUTPUT', 18X,'TIME',13X,'OUTPUT',11X,'TIME',13X,'OUTPUT',/) 22 FORMAT(8F16.7) JD=0 IF((4*JA).LT.IO) JD=IO-(4*JA) IF(JD.EQ.0) GO TO 29 DO 28 I=1,JD 28 WRITE(3,27) T(JF+I),OUTPUT(JF+I) 27 FORMAT(96X,2F16.7) 29 IRUB=1 IK=IRE IF(IK.LT.1) IK=1 C C PLOT TIME HISTORY C CALL GENPT(T,OUTPUT,200,0) DO 24 KK=1,200 T(KK)=0.0 24 OUTPUT(KK)=0.0 IF(INRNG(3,IF,4)) GO TO 93 IF(IRE.GT.6) GO TO 23 93 IF(IRE.GT.12) GO TO 23 IF(IRE.LT.1) GO TO 23 IF(INRNG(3,IF,5)) GO TO 30 23 IRUB=2 IF(JZ.EQ.1) CALL GENPT(T,OUTPUT,200,0) RETURN END C***********************************************************************C C C C SUBROUTINE TIME: C C C C SETS TIME STEP; INSTALLS LAPLACE TRANSFORM OF SPECIFIED C C FORCING FUNCTION C C***********************************************************************C C SUBROUTINE TIME(IF,AMP,PTIME,RTIME,W,OUTPUT,T,IO,NG,ING,ROOTR,ROOT 1I,IDGM1,ICODE) IMPLICIT REAL*8(N) REAL*8 ROOTR,ROOTI,GAINDG,SCHNOK(10) DIMENSION P(10),NG(10),K(10),OUT(200),OUTPUT(200),T(200),ROOTR(10) 1,ROOTI(10),SAVE(200),SZZ(100) COMMON TDEL,TMAX,SMALL,IRE C C CHECK INPUT DATA C ICODE=4 IF(IF.EQ.5.AND.W.EQ.0.0) W=.001 IF (IF.GT.6.OR.IF.LT.0) RETURN IF(IRE.GE.1) GO TO 30 WRITE(3,1000) WRITE(3,1001) IF,AMP,W,RTIME,PTIME 1000 FORMAT(/////////,10X,'PARAMETERS TRANSFERED FROM RESPON',////) 1001 FORMAT(10X,'IF=',I2,5X,'AMP=',F10.7,5X,'W=',F10.7,5X, 1'RTIME=',F10.7,5X,'PTIME=',F10.7) IF(AMP.EQ.0.0) AMP=1.0 GAINDG=1./AMP C C DETERMINE TMAX C SMALL=1.E6 DO 1 I=1,IDGM1 ABSR=DABS(ROOTR(I)) IF (ABSR.EQ.0) GO TO 1 IF (ABSR.LT.SMALL)SMALL=ABSR 1 CONTINUE IF(IF.EQ.5) GO TO 957 SMARR=0.0 JJX=0 DO 955 J=1,IDGM1 SCHNOK(J)=DBLE(ABS(ROOTI(J)))/(DBLE(ABS(ROOTR(J)))+1.0D-6) IF(SCHNOK(J).GT.SMARR) JJX=J 955 IF(SCHNOK(J).GT.SMARR) SMARR=SCHNOK(J) IF(SMARR.GT.50.0) SMALL=ABS(ROOTI(JJX))/4.0 957 TMAX=6./SMALL TMAXX=TMAX TDEL=TMAX/200 WRITE(3,1002) WRITE(3,1003) SMALL 1002 FORMAT(/,10X,'THE VALUE OF THE LONGEST TIME CONSTANT IS') 1003 FORMAT(50X,F10.7) WRITE(3,1004) WRITE(3,1005) TDEL 1004 FORMAT(10X,/,10X,'THE VALUE OF THE TIME INCREMENT IS',/) 1005 FORMAT(40X,D17.10,/) IF(IF.NE.5) GO TO 30 TMAX=12.56/W IF(TMAX.LT.5) TMAX=25.12/W 30 IF(IRE.GE.1) TMAX=TMAX+199.0*TDEL IF(IRE.EQ.0.AND.TMAX.GT.TMAXX) GO TO 400 IF((TMAX-100.*TDEL).GT.TMAXX) ICODE=0 IF((TMAX-100.*TDEL).GT.TMAXX.AND.IF.NE.5) RETURN 400 CONTINUE C C GO TO SECTION FOR SPECIFIED INPUT FORCING FUNCTION. C GO TO (2,2,4,4,14,14),IF CALL TYME(OUTPUT,T,IO,NG,ING,ROOTR,ROOTI,IDGM1,GAINDG,ICODE) RETURN C------------------------------------------------------------------------- C C THIS SECTION PREPARES THE FORCING FUNCTION FOR A C STEP INPUT. C C----------------------------------------------------------------------- 2 IDGM1=IDGM1+1 ROOTR(IDGM1)=0. ROOTI(IDGM1)=0. CALL TYME(OUTPUT,T,IO,NG,ING,ROOTR,ROOTI,IDGM1,GAINDG,ICODE) RETURN C----------------------------------------------------------------------- C 0630 C THIS SECTION PREPARES THE FORCING FUNCTION FOR A 0631 C PULSE EXCITATION OF THE SYSTEM. 0632 C 0633 C----------------------------------------------------------------------- 0634 4 IF(IRE.GE.1) GO TO 205 0635 IDGM1=IDGM1+1 0636 ROOTR(IDGM1)=0.0 0637 ROOTI(IDGM1)=0.0 0638 199 ISHO=0 0639 IF(TDEL.GE.0.025*PTIME) ISHO=1 0640 IF(ISHO.NE.1) GO TO 198 0641 TDEL=TDEL/2.0 0642 TMAX=TMAX/2.0 0643 GO TO 199 0644 198 ISHU=0 0645 IF(PTIME.GE.41.*TDEL) ISHU=1 0646 IF(ISHU.NE.1) GO TO 201 0647 TDEL=TDEL*2.0 0648 TMAX=TMAX*2.0 0649 GO TO 198 0650 201 GAINDG=GAINDG 0651 205 CALL TYME(SAVE,T,IO,NG,ING,ROOTR,ROOTI,IDGM1,GAINDG,ICODE) 0652 IF (ICODE.NE.0) GO TO 8 0653 MTIME=PTIME/TDEL+0.5 0654 IF (MTIME.EQ.0) GO TO 6 0655 DO 5 I=1,MTIME 0656 5 OUTPUT(I)=SAVE(I) 0657 6 IP1=MTIME+1 0658 DO 7 I=IP1,IO 0659 7 OUTPUT(I)=SAVE(I)-SAVE(I-MTIME) 0660 IF(IRE.EQ.0) GO TO 505 0661 DO 502 I=1,MTIME 0662 502 OUTPUT(I)=SAVE(I)-SZZ(I) 0663 505 DO 503 I=1,MTIME 0664 ICK=IO-MTIME+I-1 0665 503 SZZ(I)=SAVE(ICK) 0666 8 RETURN C----------------------------------------------------------------------- 0754 C 0755 C THIS SECTION PREPARES THE FORCING FUNCTION FOR A 0756 C SINUSOIDAL EXCITATION OF THE SYSTEM. 0757 C 0758 C----------------------------------------------------------------------- 0759 14 IF(IRE.GE.1) GO TO 333 0760 IDGM1=IDGM1+1 0761 ROOTR(IDGM1)=0. 0762 ROOTI(IDGM1)=W 0763 IDGM1=IDGM1+1 0764 ROOTR(IDGM1)=0. 0765 ROOTI(IDGM1)=-W 0766 DO 15 IG=1,ING 0767 15 NG(IG)=NG(IG)*W 0768 333 CALL TYME(OUTPUT,T,IO,NG,ING,ROOTR,ROOTI,IDGM1,GAINDG,ICODE) 0769 IF(IRE.GE.1) GO TO 334 0770 DO 16 IG=1,ING 0771 16 NG(IG)=NG(IG)/W 0772 334 RETURN END C***********************************************************************C C C C SUBROUTINE TYME C C C C FINDS POLE RESIDUES AND EVALUATES THE TIME RESPONSE C C***********************************************************************C C SUBROUTINE TYME(OUTPUT,T,IO,NG,ING,ROOTR,ROOTI,IDGM1,GAINDG,ICODE) IMPLICIT REAL*8 (C,K,L),REAL*8(N) REAL*8 P,S,OUT1,QT,QT1,QTT1,QTT2 REAL*8 ROOTR,ROOTI,GAINDG DIMENSION P(10,2),NG(10),K(10,2),OUTPUT(200),T(200),ROOTR(10), 1 ROOTI(10),TTEST(16),S(2),OUT1(2),CAINDG(2),KJ(2),QT(2), 2 QT1(2),CEXPNT(2),KNUM(2),CT(2),CUT(2) COMMON TDEL,TMAX,SMALL,IRE COMMON/FIXX/TST DATA TTEST/.001,.0025,.005,.01,.025,.05,.1,.25,.5,1.,2.5,5.,10.,25 1.,50.,100./ C C CHECK FOR BAD ENTRY C IF (IDGM1.LT.ING) GO TO 10 C C CHECK FOR MULTIPLE ROOTS C MZAP=0 0874 CAINDG(1)=GAINDG 0875 CAINDG(2)=0.0D0 0876 EPP=0.0000056 0877 1005 DO 1 I=1,IDGM1 0878 RRP1=ROOTR(I)+EPP 0879 RRM1=ROOTR(I)-EPP 0880 RIP1=ROOTI(I)+EPP 0881 RIM1=ROOTI(I)-EPP 0882 DO 1 J=1,IDGM1 0883 IF (I.EQ.J) GO TO 1 0884 RRJ=ROOTR(J) 0885 RIJ=ROOTI(J) 0886 IF (RRM1.LT.RRJ.AND.RRP1.GT.RRJ.AND.RIM1.LT.RIJ.AND.RIP1.GT.RIJ) G 0887 1O TO 9 0888 1 CONTINUE 0889 ICODE=0 0890 C C SET UP FOR SUCCESSIVE PLOT SEGMENTS C IF(IRE.GE.2) TST=TST+199.*TDEL 0894 IF(IRE.GE.1) GO TO 30 0895 TDEL=TMAX/200 0896 DO 2 I=1,16 0897 IF (TDEL.LT.TTEST(I)) GO TO 3 0898 2 CONTINUE 0899 I=16 0900 3 TDEL=TTEST(I) 0901 IF(IRE.NE.0) GO TO 11 0902 TMAX=199.*TDEL 0903 TST=TMAX 0904 11 DO 4 I=1,IDGM1 0905 P(I,1)=ROOTR(I) 0906 4 P(I,2)=ROOTI(I) 0907 C C DETERMINE THE K'S 0909 C RAT=1.0E+14 0911 DOG=0.0 0912 DO 6 J=1,IDGM1 0913 S(1)=P(J,1) 0914 S(2)=P(J,2) 0915 WRITE(3,200) S(1),S(2),NG(J) 0916 200 FORMAT(10X,'S=',D15.8,3X,D15.8,3X,'NG=',D15.8,3X) 0917 CALL CQPVAL(KNUM,S,NG,ING) 0918 KJ(1)=1.0D0 0919 KJ(2)=0.0D0 0920 DO 5 JJ=1,IDGM1 0921 IF (JJ.EQ.J) GO TO 5 0922 QTT1=-P(JJ,1) 0923 QTT2=-P(JJ,2) 0924 QT(1)=QTT1 0925 QT(2)=QTT2 0926 CALL QCADD (S,QT,QT1) 0927 CALL QCDIV (KJ,QT1,QT) 0928 QTT1=QT(1) 0929 QTT2=QT(2) 0930 KJ(1)=QTT1 0931 KJ(2)=QTT2 0932 WRITE(3,101) KJ(1),KJ(2) 0933 101 FORMAT(15X,'KJ=',D15.8,3X,D15.8) 0934 5 CONTINUE 0935 IF(CAINDG(1).NE.0.0D0) GOTO 1010 0936 CAINDG(1)=1.0D0 0937 CAINDG(2)=0.0D0 0938 1010 CONTINUE 0939 CALL QCDIV (KNUM,CAINDG,QT) 0940 CALL QCMULT (KJ,QT,QT1) 0941 K(J,1)=QT1(1) 0942 K(J,2)=QT1(2) 0943 IF(DABS(K(J,1)).LT.RAT) RAT=DABS(K(J,1)) 0944 IF(DABS(K(J,2)).GT.DOG) DOG=DABS(K(J,2)) 0945 WRITE(3,100) K(J,1),K(J,2), 0946 1 KNUM(1),KNUM(2),CAINDG(1),CAINDG(2) 0947 100 FORMAT(10X,'K(J)=',D15.8,3X,D15.8,3X,' KNUM=',D15.8,3X,D15.8, 0948 13X,'CAINDG=',D12.5,2X,D12.5) 0949 6 CONTINUE 0950 C C CHECK FOR EXCESSIVE VARIATIONS IN RESIDUE VALUES. C IF(DOG/RAT.GT.1.0E+8) WRITE(3,40) 40 FORMAT(//,30X,'WARNING: RESULTS MAY BE DISTORTED BY WIDE ', 1'DIFFERENCES IN RESIDUE VALUES. ', 2//) C C DETERMINE THE TIME RESPONSE C T1=-TDEL JCODE=0 30 IF(IRE.GE.1) T1=TST-TDEL IO=0 WRITE(3,31) IRE,T1 31 FORMAT(/,3X,'AT THE BEGINNING OF THE',I2,1X,'ITERATION, THE VALUE 1OF THE TIME IS',D14.7) WRITE(3,32) TMAX 32 FORMAT(3X,'AT THIS TIME TMAX =',D17.10,/) NBUMP=DABS(NG(1)) DO 1001 JJ=2,ING 1001 IF(DABS(NG(JJ)).GT.NBUMP) NBUMP=DABS(NG(JJ)) 7 IO=IO+1 OUT1(1)=0.0D0 OUT1(2)=0.0D0 T1=T1+TDEL DO 8 J=1,IDGM1 IF (ROOTR(J)*T1.LT.-100.) GO TO 8 LBIG=P(J,1) IF(DABS(LBIG).LT.1.0D-37)LBIG=0.0D0 LSML=P(J,2) IF(DABS(LSML).LT.1.0D-37) LSML=0.0D0 P(J,1)=LBIG P(J,2)=LSML LMOX=K(J,1) LNXT=K(J,2) C********************************************************************* C SPECIAL SECTION ADDED TO PERMIT COMPUTATION ON MACHINES C WITHOUT THE CAPABILITY OF PERFORMAING EXTENDED PRECISION C COMPLEX ARITHMETIC. C********************************************************************* CT(1)=T1 CT(2)=0.0D0 QT1(1)=P(J,1) QT1(2)=P(J,2) CALL QCMULT (CT,QT1,QT) CEXPNT(1)=DEXP(QT(1)) CEXPNT(2)=CEXPNT(1)*DSIN(QT(2)) CEXPNT(1)=CEXPNT(1)*DCOS(QT(2)) LSON=CEXPNT(1) 1000 LGRL=CEXPNT(2) 1001 C********************************************************************* 1002 C END SPECIAL SECTION 1003 C********************************************************************* 1004 IF(DABS(LSON).LT.1.0D-37) LSON=0.0D0 1005 IF(DABS(LGRL).LT.1.0D-37) LGRL=0.0D0 1006 IF(DABS(LMOX).LT.1.0D-37) LMOX=0.0D0 1007 IF(DABS(LNXT).LT.1.0D-37) LNXT=0.0D0 1008 K(J,1)=LMOX 1009 K(J,2)=LNXT 1010 CEXPNT(1)=LSON 1011 CEXPNT(2)=LGRL 1012 C********************************************************************* 1013 C THIS SECTION ADDED TO PERMIT COMPUTATION ON MACHINES 1014 C WITHOUT THE CAPABILITY OF EXTENDED PRECISION COMPLEX 1015 C ARITHMETIC 1016 C********************************************************************* 1017 CT(1)=K(J,1) 1018 CT(2)=K(J,2) 1019 CALL QCMULT (CT,CEXPNT,QT) 1020 OUT1(1)=OUT1(1)+QT(1) 1021 OUT1(2)=OUT1(2)+QT(2) 1022 C********************************************************************* 1023 C END SPECIAL SECTION 1024 C********************************************************************* 1025 8 CONTINUE 1026 ABLE=OUT1(1) 1027 OUTPUT(IO)=ABLE 1028 C C CHECK FOR EXCESSIVE COMPLEX PARTS OF OUTPUT VECTOR. C RREAL=OUT1(1) 1032 UNREAL=OUT1(2) 1033 IF(ABS(RREAL).LT.1.0E-5) GO TO 1004 1034 IF(ABS(UNREAL)/ABS(RREAL).LT.1.0E-5) GO TO 1003 1035 1004 IF (ABS(UNREAL).GT.DABS(.001D0*NBUMP)) ICODE=1 1036 IF(ICODE.EQ.1) WRITE(3,1002) UNREAL,IO 1037 IF(ICODE.EQ.1) WRITE(3,1008) RREAL 1038 1008 FORMAT(20X,'THE REAL PART OF THE OUTPUT VECTOR IS',F17.7) 1039 IF(ABS(UNREAL).LT.0.01*NBUMP) ICODE=0 1040 IF(ABS(UNREAL).GT..001*NBUMP.AND.ABS(UNREAL).LT..01*NBUMP) JCODE=1 1041 1002 FORMAT(15X,'THE COMPLEX PART OF THE OUTPUT VECTOR IS',F17.7,2X,I3) 1042 1003 CONTINUE 1043 T(IO)=T1 1044 IF(IO.EQ.200) GO TO 1006 1045 IF (T1.LT.TMAX) GO TO 7 1046 1006 IF(JCODE.EQ.1) WRITE(3,1007) 1047 1007 FORMAT(/,20X,'WARNING: OUTPUT VALUES ARE SUSPECT.',/, 1048 125X,'POSSIBLE CAUSE: REAL PARTS OF COMPLEX CONJUGATE ROOTS HAVE ' 1049 2'BEEN SEPARATED.',/) 1050 RETURN 1051 C C RELAX ICODE=2 TERMINATION CRITERION SLIGHTLY 1053 C 9 EPP=EPP-0.0000001 1055 IF(MZAP.GT.3) ICODE=2 1056 IF(MZAP.GT.3) RETURN 1057 MZAP=MZAP+1 1058 GO TO 1005 1059 10 ICODE=3 1060 RETURN 1061 END 1062 C***********************************************************************C 1063 C C 1065 C SUBROUTINE CQPVAL C 1066 C C 1067 C***********************************************************************C 1077 C SUBROUTINE CQPVAL(RES,ARG,X,IDIMX) 1079 REAL*8 RES,ARG,QT,RUS 1080 REAL*8 X 1081 DIMENSION X(20),RES(2),ARG(2),QT(2),RUS(2) 1082 C C MODIFIED TO WORK WITH QCMULT AND QCADD 1084 C RES(1)=0.0D0 1086 RES(2)=0.0D0 1087 J=IDIMX 1088 1 IF (J) 3,3,2 1089 2 QT(1)=X(J) 1090 QT(2)=0.0D0 1091 CALL QCMULT (RES,ARG,RUS) 1092 CALL QCADD (RUS,QT,RES) 1093 J=J-1 1094 GO TO 1 1095 3 RETURN 1096 END 1097 C********************************************************************** C C SUBROUTINE QCADD C C********************************************************************** SUBROUTINE QCADD(A,B,C) IMPLICIT REAL*8 (A,B,C) DIMENSION A(2),B(2),C(2) C(1)=A(1)+B(1) C(2)=A(2)+B(2) RETURN END C********************************************************************** C C SUBROUTINE QCMULT C C********************************************************************** SUBROUTINE QCMULT(A,B,C) IMPLICIT REAL*8 (A,B,C) DIMENSION A(2),B(2),C(2) C(1)=A(1)*B(1)-(A(2)*B(2)) C(2)=A(1)*B(2)+(A(2)*B(1)) RETURN END C********************************************************************** C C SUBROUTINE QCDIV C C********************************************************************** SUBROUTINE QCDIV(A,B,C) IMPLICIT REAL*8 (A,B,C) DIMENSION A(2),B(2),C(2) QD=B(1)**2+B(2)**2 C(1)=(A(1)*B(1)+A(2)*B(2))/QD C(2)=(A(2)*B(1)-A(1)*B(2))/QD RETURN END C********************************************************************** C C SUBROUTINE LOCUS C C PROVIDES CONTROL FOR THE GENERATION OF ROOT LOCUS C PLOTS WHEN REQUESTED; OTHERWISE, FINDS ROOTS OF C CHARACTERISTIC EQUATION C********************************************************************** SUBROUTINE LOCUS(A,AK,Z,DK,N,M,MM1,MM2,JK,ISW) 2751 IMPLICIT REAL*8(A-H,O-Z) 2752 DIMENSION A(11),AK(11),BZ(11),ZF(20),AC(11),B(11),Z(20) 2753 DIMENSION ZZ(20),ZEROR(10),ZEROI(10),TH(10),SIGN(10) 2754 COMMON/LOCO/RR(200,10),RI(200,10) 2755 C 2756 MM1=M-1 2757 MM2=2*MM1 2758 IG=0 2759 IF(N.EQ.-1) GO TO 31 2760 32 JK=0 2761 DO 77 KK=1,10 2762 77 AC(KK)=0.0D0 2763 C 2764 C ZERO OUT RR AND RI ARRAYS. 2765 C 2766 DO 3 LL=1,N 2767 DO 3 LM=1,9 2768 RR(LL,LM)=0.0D0 2769 3 RI(LL,LM)=0.0D0 2770 DO 4 J=1,N 2771 BK=DK*DBLE(FLOAT(J-1)) 2772 DO 5 I=1,M 2773 AC(I)=A(I)+AK(I)*DK*DBLE(FLOAT(J-1)) 2774 5 CONTINUE 2775 C 2776 C REVERSE ORDER OF COEFFICIENTS FOR INPUT TO ZRPOLY. 2777 C 2778 DO 13 I=1,M 2779 MP=M+1-I 2780 13 B(MP)=AC(I) 2781 MP1=M+1 2782 DO 25 I=MP1,11 2783 25 B(I)=0.0D0 2784 C 2785 IF(ISW.NE.4) GO TO 47 2786 IF(ISW.EQ.4) ISW=0 2787 WRITE(3,44) 2788 44 FORMAT(/,5X,'THE INPUTS TO ZRPOLY ARE:',/) 2789 WRITE(3,45) (B(LMM),LMM=1,M) 2790 45 FORMAT(2X,5(D23.16,2X)) 2791 47 CONTINUE 2792 C 2793 C SWITCH TO FORMULATION FOR PURE IMAGINARY ROOTS 2794 C 2795 IZK=0 2796 DO 15 JX=2,10,2 2797 15 IF(B(JX).EQ.0.0D0) IZK=IZK+1 2798 IF(IZK.NE.5) GO TO 16 2799 WRITE(3,24) 2800 24 FORMAT(103X,'ALL ROOTS ARE IMAGINARY') 2801 JC=0 2802 DO 17 JX=1,11,2 2803 JC=JC+1 2804 17 BZ(JC)=B(JX) 2805 MX=MP1/2 2806 MX1=MX-1 2807 WRITE(3,26) 2808 26 FORMAT(/,30X,'INPUTS FOR THE REDUCED FORM OF THE CHARACT', 2809 1'ERISTIC EQUATION ARE:',/) 2810 WRITE(3,27) (BZ(MQ),MQ=1,MX) 2811 27 FORMAT(1X,40X,F20.12) 2812 C 2813 CALL ZRPOLY(BZ,MX1,ZF,IER) 2814 C 2815 WRITE(3,28) 2816 28 FORMAT(/,30X,'ROOTS OF THIS REDUCED EQUATION ARE:',/) 2817 MX1T2=MX1*2 2818 MX1T1=MX1T2-1 2819 DO 35 MQ=1,MX1T1,2 2820 DIFF=ZF(MQ)-ZF(MQ+2) 2821 35 IF(DBLE(ABS(DIFF/ZF(MQ))).LT.1.0D-7.AND.DBLE(ABS(DIFF)).NE.0.0D0) 2822 1CALL NEWTON(BZ,ZF,MQ,MX) 2823 WRITE(3,29) (ZF(MQ),MQ=1,MX1T2) 2824 29 FORMAT(1X,40X,F20.12,' + J',F20.12) 2825 DO 20 JX=1,19,2 2826 20 Z(JX)=0.0D0 2827 DO 21 I=1,MX1T2,2 2828 LOM=2*I+2 2829 K=2*I 2830 Z(K)=DSQRT(DBLE(ABS(ZF(I)))) 2831 21 Z(LOM)=-DSQRT(DBLE(ABS(ZF(I)))) 2832 LOM=4*MX1+1 2833 DO 22 JX=LOM,20 2834 22 Z(JX)=0.0D0 2835 GO TO 23 2836 C 2837 C CALL ROOT SOLVING ROUTINE 2838 C 2839 16 CALL ZRPOLY(B,MM1,Z,IER) 2840 C 2841 DO 14 JJ=1,MM2 2842 14 IF(DBLE(ABS(Z(JJ))).LT.1.0D-15) Z(JJ)=0.0D0 2843 23 IF(ISW.EQ.1) GO TO 1 2844 C 2845 C WRITE OUT RESULTS: 2846 C INCLUDE DAMPING RATIO AND NATURAL FREQUENCY 2847 C IN SECOND GROUP OF EACH PRINTOUT 2848 C 2849 WRITE(3,6) BK 2850 WRITE(3,7) (Z(JJ),JJ=1,MM2) 2851 1 CONTINUE 2852 JK=JK+1 2853 JN=0 2854 WRITE(3,30) 2855 30 FORMAT(/,48X,'ROOT',6X,'DAMPING RATIO',6X,'NATURAL FREQUENCY') 2856 DO 8 JL=1,MM2,2 2857 JN=JN+1 2858 IF(Z(JL).EQ.0.0D0.AND.Z(JL+1).EQ.0.0D0) GO TO 18 2859 IF(Z(JL).GT.0.0D0) GO TO 18 2860 ZETA=DBLE(ABS(Z(JL)))/DSQRT(Z(JL)*Z(JL)+Z(JL+1)*Z(JL+1)) 2861 IF(ZETA.EQ.0.0D0) OMEGA=Z(JL+1) 2862 IF(ZETA.EQ.0.0D0) GO TO 19 2863 OMEGA=DBLE(ABS(Z(JL)))/ZETA 2864 IF(ZETA.EQ.1.0D0) OMEGA=0.0D0 2865 19 WRITE(3,2) JN,ZETA,OMEGA 2866 2 FORMAT(47X,I4,3X,F17.10,3X,F17.10) 2867 18 CONTINUE 2868 RR(JK,JN)=Z(JL) 2869 8 RI(JK,JN)=Z(JL+1) 2870 IF(J.GT.1) GO TO 4 2871 POLE=0.0D0 2872 DO 106 JJ=1,MM2,2 2873 106 POLE=POLE+Z(JJ) 2874 4 CONTINUE 2875 IF(N.LT.2) RETURN 2876 WRITE(3,12) 2877 12 FORMAT(1H1) 2878 C 2879 C CALL PLOTTING ROUTINE 2880 C 2881 CALL PPLOT(N,MM1) 2882 DO 100 J=1,11 2883 100 B(J)=0.0D0 2884 DO 101 J=1,M 2885 MP=M+1-J 2886 101 IF(AK(MP).NE.0.0D0) GO TO 102 2887 102 CONTINUE 2888 C 2889 C DETERMINE ROOTS OF NUMERATOR 2890 C 2891 DO 103 J=1,MP 2892 MR=MP+1-J 2893 103 B(MR)=AK(J) 2894 MMR=MP-1 2895 IF(MP.EQ.1) RETURN 2896 DO 111 J=1,10 2897 MPP=11-J 2898 111 IF(AK(MPP).NE.0.0D0) GO TO 112 2899 112 IF(MPP.GT.M) WRITE(3,113) 2900 113 FORMAT(//,5X,'NUMBER OF TERMS IN H-FUNCTION GREATER', 2901 1' THAN IDG.',/,5X,'PROBABLE USER ERROR',//) 2902 IF(MPP.GT.M) RETURN 2903 IF(MMR.EQ.MM1) WRITE(3,314) 2904 314 FORMAT(//,5X,'NUMBER OF ZEROS EQUAL NUMBER OF POLES',/, 2905 15X,'NO ASYMPTOTES',//) 2906 IF(MMR.EQ.MM1) RETURN 2907 MMR2=MMR*2 2908 CALL ZRPOLY(B,MMR,ZZ,IER) 2909 JN=0 2910 DO 104 J=1,MMR2,2 2911 JN=JN+1 2912 ZEROR(JN)=ZZ(J) 2913 104 ZEROI(JN)=ZZ(J+1) 2914 WRITE(3,117) 2915 117 FORMAT(1H1,////////,15X,'THE ROOTS OF THE NUMERATOR ARE:',/) 2916 WRITE(3,118) (ZEROR(JMB),ZEROI(JMB),JMB=1,MMR) 2917 118 FORMAT(15X,F23.15,2X,'+J',F23.15) 2918 WRITE(3,131) 2919 131 FORMAT(/,15X,'NUMERATOR ROOTS WILL BE CORRECT ONLY', 2920 1/,15X,'FOR UNITY FEEDBACK. FOR OTHER "H" FUNCTIONS',/, 2921 215X,'THESE RESULTS SHOULD BE IGNORED...FOR UNITY ',/, 2922 315X,'FEEDBACK ALSO, THE FREQUENCY AT WHICH THE BODE',/, 2923 415X,'PLOT HAS A ZERO PHASE MARGIN CAN ALSO BE USED TO',/, 2924 515X,'DETERMINE THE GAIN MARGIN OR VALUE OF THE FEEDBACK', 2925 6/,15X,'GAIN AT WHICH THE ROOT LOCUS WILL CROSS THE ',/, 2926 715X,'IMAGINARY AXIS, I.E., GO UNSTABLE.') 2927 ZERO=0.0D0 2928 C 2929 C DETERMINE CG OF ASYMPTOTES 2930 C 2931 DO 105 J=1,MMR 2932 105 ZERO=ZERO+ZEROR(J) 2933 PI=DACOS(-1.0D0) 2934 RCG=(POLE-ZERO)/DBLE(FLOAT(MM1-MMR)) 2935 C 2936 C DETERMINE ASYMPTOTE ANGLES 2937 C 2938 K=MM1-MMR 2939 IF(K.LE.0) RETURN 2940 DO 107 J=1,K 2941 TH(J)=(2.0D0*DBLE(FLOAT(J))-1.0D0)/DBLE(FLOAT(MM1-MMR))*PI 2942 SIGN(J)=1.0D0 2943 107 IF(TH(J).GT.(PI/2.0D0).AND.TH(J).LT.(1.5D0*PI)) SIGN(J)=-1.0D0 2944 WRITE(3,119) POLE,ZERO 2945 119 FORMAT(/,5X,'SUM OF POLES=',F10.3,2X,'SUM OF ZEROS=',F10.3,/) 2946 WRITE(3,115) RCG 2947 115 FORMAT(//,15X,'ROOT CENTER OF GRAVITY=',F23.15,/) 2948 WRITE(3,116) (TH(J),J=1,K) 2949 116 FORMAT(15X,'ANGLE IN RADIANS OF ASYMPTOTE=',F15.10) 2950 C 2951 C GENERATE ASYMPTOTES 2952 C 2953 DO 108 I=1,200 2954 DO 108 J=1,10 2955 RR(I,J)=0.0D0 2956 108 RI(I,J)=0.0D0 2957 DO 109 I=1,K 2958 DO 109 J=1,100 2959 RR(J,I)=RCG+DBLE(ABS(RCG))/50.0D0*DBLE(FLOAT(J))*SIGN(I) 2960 IF(TH(I).EQ.PI/2.0D0.OR.TH(I).EQ.1.5D0*PI) RR(J,I)=RCG 2961 SIGN1=1.0D0 2962 IF(TH(I).GT.PI/2.0D0.AND.TH(I).LT.1.5D0*PI) SIGN1=-1.0D0 2963 IF(TH(I).EQ.PI/2.0D0.OR.TH(I).EQ.1.5D0*PI) TH(I)=TH(I)/(1.0D0 2964 1+1.0D-5) 2965 DONT=DTAN(TH(I))*SIGN1 2966 IF(DONT.GT.1.0D+4) DONT=1.0D+4 2967 IF(DONT.LT.-1.0D+4) DONT=-1.0D+4 2968 109 RI(J,I)=DBLE(ABS(RCG))/50.0D0*DBLE(FLOAT(J))*DONT 2969 WRITE(3,12) 2970 CALL PPLOT(100,K) 2971 WRITE(3,110) 2972 110 FORMAT(15X,'PLOTTED ABOVE ARE THE ASYMPTOTES OF THE LOCUS') 2973 6 FORMAT(//,10X,'THE VALUE OF THE FEEDBACK GAIN IS',2X,D17.10,/) 2974 7 FORMAT(30X,D17.10,3X,'+J',D17.10) 2975 IF(IG.GT.0) GO TO 31 2976 RETURN 2977 C 2978 C IF N=-1, SETUP FOR AUTOMATIC SELECTION OF ROOT 2979 C LOCUS INTERVAL AND NUMBER OF POINTS ON LOCUS 2980 C 2981 31 IG=IG+1 2982 IF(IG.EQ.3) RETURN 2983 DK=0.05D0*(A(1)+AK(1))**(1.0D0/DBLE(FLOAT(M-1))) 2984 N=40 2985 IF(IG.EQ.2) DK=10.0D0*DK 2986 WRITE(3,34) 2987 34 FORMAT(/,25X,'AUTOMATIC INCREMENT SELECTION ALGORITHM ', 1'INVOKED') WRITE(3,33) A(1),AK(1),M,DK 33 FORMAT(//,10X,'NON-CHANGING CONSTANT IN EQUATION=',F17.10,/, 110X,'CHANGING CONSTANT IN EQUATION=',F17.10,/, 210X,'NUMBER OF TERMS IN CHARACTERISTIC EQUATION=',I4,/, 310X,'INCREMENT IN GAIN SELECTED BY PROGRAM=',F17.10,/) GO TO 32 END C********************************************************************** C C SUBROUTINE ZRPOLY C C MODIFICATION OF 1978 IMSL ROUTINE TO FIND THE ROOTS OF C POLYNOMIALS C********************************************************************** SUBROUTINE ZRPOLY (A,NDEG,Z,IER) INTEGER NDEG,IER DOUBLE PRECISION A(11),Z(20) INTEGER N,NN,J,JJ,I,NM1,ICNT,N2,L,NZ,NPI REAL*8 ETA,RMRE,RINFP,REPSP,RADIX,RLO,XX,YY,SINR, 1 COSR,RMAX,RMIN,X,SC,XM,FF,DX,DF,BND,XXX,ARE REAL*8 PT(101),RASH DOUBLE PRECISION TEMP(101),P(101),QP(101),RK(101),QK(101), 1 SVK(101) DOUBLE PRECISION SR,SI,U,V,RA,RB,C,D,A1,A2,A3, 1 A6,A7,E,F,G,H,SZR,SZI,RLZR,RLZI, 2 T,AA,BB,CC,FACTOR,REPSR1,ZERO,ONE,FN LOGICAL ZEROK CHARACTER *6 NAME COMMON /ZRPQLJ/ P,QP,RK,QK,SVK,SR,SI,U,V,RA,RB,C,D,A1,A2,A3,A6, 1 A7,E,F,G,H,SZR,SZI,RLZR,RLZI,ETA,ARE,RMRE,N,NN DATA NAME/'ZRPOLY'/ DATA RINFP/1.000D+60/ DATA REPSP/1.00D-60/ DATA RADIX/16.0D0/ DATA REPSR1/1.000000000000D-16/ DATA ZERO/0.0D0/,ONE/1.0D0/ C IER = 0 IF (NDEG .GT. 100 .OR. NDEG .LT. 1) GO TO 165 ETA = REPSR1 ARE = ETA RMRE = ETA RLO = REPSP/ETA XX=0.5D0*DSQRT(2.0D0) RASH=94.0D0/180.0D0*DACOS(-1.00D0) SINR=DSIN(RASH) COSR=DCOS(RASH) YY = -XX N = NDEG NN = N+1 IF (A(1).NE.ZERO) GO TO 5 IER = 130 GO TO 9000 5 IF (A(NN).NE.ZERO) GO TO 10 J = NDEG-N+1 JJ = J+NDEG Z(J) = ZERO Z(JJ) = ZERO NN = NN-1 N = N-1 IF (NN.EQ.1) GO TO 9005 GO TO 5 10 DO 15 I=1,NN P(I) = A(I) 15 CONTINUE 20 IF (N.GT.2) GO TO 30 IF (N.LT.1) GO TO 9005 IF (N.EQ.2) GO TO 25 Z(NDEG) = -P(2)/P(1) Z(NDEG+NDEG) = ZERO GO TO 145 25 CALL ZRPQLI (P(1),P(2),P(3),Z(NDEG-1),Z(NDEG+NDEG-1),Z(NDEG), 1 Z(NDEG+NDEG)) GO TO 145 30 RMAX =0.D0 RMIN = RINFP DO 35 I=1,NN X = DBLE(ABS( (P(I)))) IF (X.GT.RMAX) RMAX = X IF (X.NE.0.D0.AND.X.LT.RMIN) RMIN = X 35 CONTINUE SC = RLO/RMIN IF (SC.GT.1.D0) GO TO 40 IF (RMAX.LT.10.D0) GO TO 55 IF (SC.EQ.0.D0) SC=REPSP*RADIX*RADIX GO TO 45 40 IF (RINFP/SC.LT.RMAX) GO TO 55 45 L = DLOG(SC)/DLOG(RADIX)+.5 IF (L .EQ. 0) GO TO 55 FACTOR = RADIX**L DO 50 I=1,NN 50 P(I) = FACTOR*P(I) 55 DO 60 I=1,NN 60 PT(I) = DBLE(ABS( (P(I)))) PT(NN) = -PT(NN) X =DEXP((DLOG(-PT(NN))-DLOG(PT(1)))/N) IF (PT(N).EQ.0.D0) GO TO 65 XM = -PT(NN)/PT(N) IF (XM.LT.X) X = XM 65 XM = X*.1 FF = PT(1) DO 70 I=2,NN 70 FF = FF*XM+PT(I) IF (FF.LE.0.D0) GO TO 75 X = XM GO TO 65 75 DX = X 80 IF (DBLE(ABS(DX/X)).LE..005D0) GO TO 90 FF = PT(1) ZRPA3200 DF = FF ZRPA3201 DO 85 I=2,N ZRPA3202 FF = FF*X+PT(I) ZRPA3203 DF = DF*X+FF ZRPA3204 85 CONTINUE ZRPA3205 FF = FF*X+PT(NN) ZRPA3206 DX = FF/DF ZRPA3207 X = X-DX ZRPA3208 GO TO 80 ZRPA3209 90 BND = X ZRPA3210 NM1 = N-1 ZRPA3214 FN = ONE/N ZRPA3215 DO 95 I=2,N ZRPA3216 95 RK(I) = (NN-I)*P(I)*FN ZRPA3217 RK(1) = P(1) ZRPA3218 AA = P(NN) ZRPA3219 BB = P(N) ZRPA3220 ZEROK = RK(N).EQ.ZERO ZRPA3221 DO 115 JJ=1,5 ZRPA3222 CC = RK(N) ZRPA3223 IF (ZEROK) GO TO 105 ZRPA3224 T = -AA/CC ZRPA3227 DO 100 I=1,NM1 ZRPA3228 J = NN-I ZRPA3229 RK(J) = T*RK(J-1)+P(J) ZRPA3230 100 CONTINUE ZRPA3231 RK(1) = P(1) ZRPA3232 ZEROK = DBLE(ABS(RK(N))).LE.DABS(BB)*ETA*10.D0 ZRPA3233 GO TO 115 ZRPA3234 105 DO 110 I=1,NM1 ZRPA3236 J = NN-I ZRPA3237 RK(J) = RK(J-1) ZRPA3238 110 CONTINUE ZRPA3239 RK(1) = ZERO ZRPA3240 ZEROK = RK(N).EQ.ZERO ZRPA3241 115 CONTINUE ZRPA3242 DO 120 I=1,N ZRPA3244 120 TEMP(I) = RK(I) ZRPA3245 DO 140 ICNT=1,20 ZRPA3248 XXX = COSR*XX-SINR*YY ZRPA3255 YY = SINR*XX+COSR*YY ZRPA3256 XX = XXX ZRPA3257 SR = BND*XX ZRPA3258 SI = BND*YY ZRPA3259 U = -SR-SR ZRPA3260 V = BND*BND ZRPA3261 CALL ZRPQLB (20*ICNT,NZ) ZRPA3264 IF (NZ.EQ.0) GO TO 130 ZRPA3265 J = NDEG-N+1 ZRPA3272 JJ = J+NDEG ZRPA3273 Z(J) = SZR ZRPA3274 Z(JJ) = SZI ZRPA3275 NN = NN-NZ ZRPA3276 N = NN-1 ZRPA3277 DO 125 I=1,NN ZRPA3278 125 P(I) = QP(I) ZRPA3279 IF (NZ.EQ.1) GO TO 20 ZRPA3280 Z(J+1) = RLZR ZRPA3281 Z(JJ+1) = RLZI ZRPA3282 GO TO 20 ZRPA3283 130 DO 135 I=1,N ZRPA3287 135 RK(I) = TEMP(I) ZRPA3288 140 CONTINUE ZRPA3289 IER = 131 ZRPA3292 145 DO 150 I=1,NDEG ZRPA3294 NPI= NDEG+I ZRPA3295 P(I) = Z(NPI) ZRPA3296 150 CONTINUE ZRPA3297 N2 = NDEG+NDEG ZRPA3298 J = NDEG ZRPA3299 DO 155 I=1,NDEG ZRPA3300 Z(N2-1) = Z(J) ZRPA3301 Z(N2) = P(J) ZRPA3302 N2 = N2-2 ZRPA3303 J = J-1 ZRPA3304 155 CONTINUE ZRPA3305 IF (IER .EQ. 0) GO TO 9005 ZRPA3306 N2 = 2*(NDEG-NN)+3 ZRPA3308 DO 160 I=1,N ZRPA3309 Z(N2) = RINFP ZRPA3310 Z(N2+1) = RINFP ZRPA3311 N2 = N2+2 ZRPA3312 160 CONTINUE ZRPA3313 GO TO 9000 ZRPA3314 165 IER = 129 ZRPA3315 9000 CONTINUE ZRPA3316 CALL UERTST (IER,NAME) ZRPA3317 9005 RETURN ZRPA3318 END ZRPA3319 SUBROUTINE ZRPQLB (L2,NZ) ZRPB3349 INTEGER L2,NZ ZRPB3351 INTEGER N,NN,J,ITYPE,I,IFLAG ZRPB3353 REAL*8 ARE,BETAS,BETAV,ETA,OSS,OTS,OTV,OVV,RMRE,SS, ZRPB3354 1 TS,TSS,TV,TVV,VV ZRPB3355 DOUBLE PRECISION P(101),QP(101),RK(101),QK(101),SVK(101) ZRPB3356 DOUBLE PRECISION SR,SI,U,V,RA,RB,C,D,A1,A2,A3, ZRPB3357 1 A6,A7,E,F,G,H,SZR,SZI,RLZR,RLZI, ZRPB3358 2 SVU,SVV,UI,VI,S,ZERO ZRPB3359 LOGICAL VPASS,SPASS,VTRY,STRY ZRPB3360 COMMON /ZRPQLJ/ P,QP,RK,QK,SVK,SR,SI,U,V,RA,RB,C,D,A1,A2,A3,A6,ZRPB3361 1 A7,E,F,G,H,SZR,SZI,RLZR,RLZI,ETA,ARE,RMRE,N,NN ZRPB3362 DATA ZERO/0.0D0/ ZRPB3363 NZ = 0 ZRPB3365 BETAV = .25D0 ZRPB3375 BETAS = .25D0 ZRPB3376 OSS = SR ZRPB3377 OVV = V ZRPB3378 CALL ZRPQLH (NN,U,V,P,QP,RA,RB) ZRPB3381 CALL ZRPQLE (ITYPE) ZRPB3382 DO 40 J=1,L2 ZRPB3383 CALL ZRPQLF (ITYPE) ZRPB3386 CALL ZRPQLE (ITYPE) ZRPB3387 CALL ZRPQLG (ITYPE,UI,VI) ZRPB3388 VV = VI ZRPB3389 SS = 0.D0 ZRPB3391 IF (RK(N).NE.ZERO) SS = -P(NN)/RK(N) ZRPB3392 TV = 1.D0 ZRPB3393 TS = 1.D0 ZRPB3394 IF (J.EQ.1.OR.ITYPE.EQ.3) GO TO 35 ZRPB3395 IF (VV.NE.0.) TV = DABS((VV-OVV)/VV) ZRP3397 IF (SS.NE.0.) TS = DABS((SS-OSS)/SS) ZRP3398 TVV = 1.D0 ZRPB3401 IF (TV.LT.OTV) TVV = TV*OTV ZRPB3402 TSS = 1.D0 ZRPB3403 IF (TS.LT.OTS) TSS = TS*OTS ZRPB3404 VPASS = TVV.LT.BETAV ZRPB3406 SPASS = TSS.LT.BETAS ZRPB3407 IF (.NOT.(SPASS.OR.VPASS)) GO TO 35 ZRPB3408 SVU = U ZRPB3412 SVV = V ZRPB3413 DO 5 I=1,N ZRPB3414 5 SVK(I) = RK(I) ZRPB3415 S = SS ZRPB3416 VTRY = .FALSE. ZRPB3419 STRY = .FALSE. ZRPB3420 IF (SPASS.AND.((.NOT.VPASS).OR.TSS.LT.TVV)) GO TO 20 ZRPB3421 10 CALL ZRPQLC (UI,VI,NZ) ZRPB3422 IF (NZ.GT.0) RETURN ZRPB3423 VTRY = .TRUE. ZRPB3428 BETAV = BETAV*.25D0 ZRPB3429 IF (STRY.OR.(.NOT.SPASS)) GO TO 25 ZRPB3433 DO 15 I=1,N ZRPB3434 15 RK(I) = SVK(I) ZRPB3435 20 CALL ZRPQLD (S,NZ,IFLAG) ZRPB3436 IF (NZ.GT.0) RETURN ZRPB3437 STRY = .TRUE. ZRPB3441 BETAS = BETAS*.25D0 ZRPB3442 IF (IFLAG.EQ.0) GO TO 25 ZRPB3443 UI = -(S+S) ZRPB3447 VI = S*S ZRPB3448 GO TO 10 ZRPB3449 25 U = SVU ZRPB3451 V = SVV ZRPB3452 DO 30 I=1,N ZRPB3453 30 RK(I) = SVK(I) ZRPB3454 IF (VPASS.AND.(.NOT.VTRY)) GO TO 10 ZRPB3458 CALL ZRPQLH (NN,U,V,P,QP,RA,RB) ZRPB3461 CALL ZRPQLE (ITYPE) ZRPB3462 35 OVV = VV ZRPB3463 OSS = SS ZRPB3464 OTV = TV ZRPB3465 OTS = TS ZRPB3466 40 CONTINUE ZRPB3467 RETURN ZRPB3468 END ZRPB3469 SUBROUTINE ZRPQLC (UU,VV,NZ) ZRPC3498 INTEGER NZ ZRPC3500 DOUBLE PRECISION UU,VV ZRPC3501 INTEGER N,NN,J,I,ITYPE ZRPC3503 REAL*8 ARE,EE,ETA,OMP,RELSTP,RMP,RMRE,T,ZM ZRPC3504 DOUBLE PRECISION P(101),QP(101),RK(101),QK(101),SVK(101) ZRPC3505 DOUBLE PRECISION SR,SI,U,V,RA,RB,C,D,A1,A2,A3, ZRPC3506 1 A6,A7,E,F,G,H,SZR,SZI,RLZR,RLZI, ZRPC3507 2 UI,VI,ZERO,PT01,ONE ZRPC3508 LOGICAL TRIED ZRPC3509 COMMON /ZRPQLJ/ P,QP,RK,QK,SVK,SR,SI,U,V,RA,RB,C,D,A1,A2,A3,A6,ZRPC3510 1 A7,E,F,G,H,SZR,SZI,RLZR,RLZI,ETA,ARE,RMRE,N,NN ZRPC3511 DATA ZERO,PT01,ONE/0.0D0,0.01D0,1.0D0/ ZRPC3512 NZ = 0 ZRPC3514 TRIED = .FALSE. ZRPC3522 U = UU ZRPC3523 V = VV ZRPC3524 J = 0 ZRPC3525 5 CALL ZRPQLI (ONE,U,V,SZR,SZI,RLZR,RLZI) ZRPC3527 IF ( DABS(DABS(SZR)-DABS(RLZR)).GT.PT01*DABS(RLZR)) RETURN ZRPC3531 CALL ZRPQLH (NN,U,V,P,QP,RA,RB) ZRPC3534 RMP = DABS(RA-SZR*RB)+DABS(SZI*RB) ZR 3535 ZM =DSQRT(DABS( (V))) ZRPC3538 EE = 2.00D0*DABS((QP(1))) ZRPC3539 T = -SZR*RB ZRPC3540 DO 10 I=2,N ZRPC3541 10 EE = EE*ZM+DABS( (QP(I))) ZRPC3542 EE = EE*ZM+DABS(SNGL(RA)+T) ZRPC3543 EE = (5.D0*RMRE+4.D0*ARE)*EE-(5.D0*RMRE+2.D0*ARE)*(DABS( (RA) ZRPC3544 1 +T)+DABS( (RB))*ZM)+2.D0*ARE*DABS(T) ZRPC3545 IF (RMP.GT.20.D0*EE) GO TO 15 ZRPC3549 NZ = 2 ZRPC3550 RETURN ZRPC3551 15 J = J+1 ZRPC3552 IF (J.GT.20) RETURN ZRPC3554 IF (J.LT.2) GO TO 25 ZRPC3555 IF (RELSTP.GT..01D0.OR.RMP.LT.OMP.OR.TRIED) GO TO 25 ZRPC3556 IF (RELSTP.LT.ETA) RELSTP = ETA ZRPC3561 RELSTP =DSQRT(RELSTP) ZRPC3562 U = U-U*RELSTP ZRPC3563 V = V+V*RELSTP ZRPC3564 CALL ZRPQLH (NN,U,V,P,QP,RA,RB) ZRPC3565 DO 20 I=1,5 ZRPC3566 CALL ZRPQLE (ITYPE) ZRPC3567 CALL ZRPQLF (ITYPE) ZRPC3568 20 CONTINUE ZRPC3569 TRIED = .TRUE. ZRPC3570 J = 0 ZRPC3571 25 OMP = RMP ZRPC3572 CALL ZRPQLE (ITYPE) ZRPC3575 CALL ZRPQLF (ITYPE) ZRPC3576 CALL ZRPQLE (ITYPE) ZRPC3577 CALL ZRPQLG (ITYPE,UI,VI) ZRPC3578 IF (VI.EQ.ZERO) RETURN ZRPC3581 RELSTP = DABS((VI-V)/VI) ZRP 3582 U = UI ZRPC3583 V = VI ZRPC3584 GO TO 5 ZRPC3585 END ZRPC3586 SUBROUTINE ZRPQLD (SSS,NZ,IFLAG) ZRPD3615 INTEGER NZ,IFLAG ZRPD3617 DOUBLE PRECISION SSS ZRPD3618 INTEGER N,NN,J,I ZRPD3620 REAL*8 ARE,EE,ETA,OMP,RMP,RMS,RMRE ZRPD3621 DOUBLE PRECISION P(101),QP(101),RK(101),QK(101),SVK(101) ZRPD3622 DOUBLE PRECISION SR,SI,U,V,RA,RB,C,D,A1,A2,A3, ZRPD3623 1 A6,A7,E,F,G,H,SZR,SZI,RLZR,RLZI, ZRPD3624 2 PV,RKV,T,S,ZERO,PT001 ZRPD3625 COMMON /ZRPQLJ/ P,QP,RK,QK,SVK,SR,SI,U,V,RA,RB,C,D,A1,A2,A3,A6,ZRPD3626 1 A7,E,F,G,H,SZR,SZI,RLZR,RLZI,ETA,ARE,RMRE,N,NN ZRPD3627 DATA ZERO/0.0D0/,PT001/0.001D0/ ZRPD3628 NZ = 0 ZRPD3636 S = SSS ZRPD3637 IFLAG = 0 ZRPD3638 J = 0 ZRPD3639 5 PV = P(1) ZRPD3641 QP(1) = PV ZRPD3643 DO 10 I=2,NN ZRPD3644 PV = PV*S+P(I) ZRPD3645 QP(I) = PV ZRPD3646 10 CONTINUE ZRPD3647 RMP = DABS(PV) ZRP 3648 RMS = DABS(S) ZRP 3651 EE = (RMRE/(ARE+RMRE))*DABS( (QP(1))) ZRP3652 DO 15 I=2,NN ZRPD3653 15 EE = EE*RMS+DABS( (QP(I))) ZRP3654 IF (RMP.GT.20.D0*((ARE+RMRE)*EE-RMRE*RMP)) GO TO 20 ZRPD3658 NZ = 1 ZRPD3659 SZR = S ZRPD3660 SZI = ZERO ZRPD3661 RETURN ZRPD3662 20 J = J+1 ZRPD3663 IF (J.GT.10) RETURN ZRPD3665 IF (J.LT.2) GO TO 25 ZRPD3666 IF (DABS(T).GT.PT001*DABS(S-T).OR.RMP.LE.OMP) GO TO 25 ZR 3667 IFLAG = 1 ZRPD3672 SSS = S ZRPD3673 RETURN ZRPD3674 25 OMP = RMP ZRPD3677 RKV = RK(1) ZRPD3680 QK(1) = RKV ZRPD3681 DO 30 I=2,N ZRPD3682 RKV = RKV*S+RK(I) ZRPD3683 QK(I) = RKV ZRPD3684 30 CONTINUE ZRPD3685 IF (DABS(RKV).LE.DABS(RK(N))*10.0D0*ETA) GO TO 40 ZR 3686 T = -PV/RKV ZRPD3690 RK(1) = QP(1) ZRPD3691 DO 35 I=2,N ZRPD3692 35 RK(I) = T*QK(I-1)+QP(I) ZRPD3693 GO TO 50 ZRPD3694 40 RK(1) = ZERO ZRPD3696 DO 45 I=2,N ZRPD3697 45 RK(I) = QK(I-1) ZRPD3698 50 RKV = RK(1) ZRPD3699 DO 55 I=2,N ZRPD3700 55 RKV = RKV*S+RK(I) ZRPD3701 T = ZERO ZRPD3702 IF (DABS(RKV).GT.DABS(RK(N))*10.0D0*ETA) T=-PV/RKV ZR 3703 S = S+T ZRPD3704 GO TO 5 ZRPD3705 END ZRPD3706 SUBROUTINE ZRPQLE (ITYPE) ZRPE3735 INTEGER ITYPE ZRPE3737 INTEGER N,NN ZRPE3739 REAL*8 ARE,ETA,RMRE ZRPE3740 DOUBLE PRECISION P(101),QP(101),RK(101),QK(101),SVK(101) ZRPE3741 DOUBLE PRECISION SR,SI,U,V,RA,RB,C,D,A1,A2,A3, ZRPE3742 1 A6,A7,E,F,G,H,SZR,SZI,RLZR,RLZI ZRPE3743 COMMON /ZRPQLJ/ P,QP,RK,QK,SVK,SR,SI,U,V,RA,RB,C,D,A1,A2,A3,A6,ZRPE3744 1 A7,E,F,G,H,SZR,SZI,RLZR,RLZI,ETA,ARE,RMRE,N,NN ZRPE3745 CALL ZRPQLH (N,U,V,RK,QK,C,D) ZRPE3757 IF (DABS(C).GT.DABS(RK(N))*100.D0*ETA) GO TO 5 ZRPE3758 IF (DABS(D).GT.DABS(RK(N-1))*100.D0*ETA) GO TO 5 ZRPE3759 ITYPE = 3 ZRPE3760 RETURN ZRPE3763 5 IF (DABS(D).LT.DABS(C)) GO TO 10 ZRPE3764 ITYPE = 2 ZRPE3765 E = RA/D ZRPE3768 F = C/D ZRPE3769 G = U*RB ZRPE3770 H = V*RB ZRPE3771 A3 = (RA+G)*E+H*(RB/D) ZRPE3772 A1 = RB*F-RA ZRPE3773 A7 = (F+U)*RA+H ZRPE3774 RETURN ZRPE3775 10 ITYPE = 1 ZRPE3776 E = RA/C ZRPE3779 F = D/C ZRPE3780 G = U*E ZRPE3781 H = V*RB ZRPE3782 A3 = RA*E+(H/C+G)*RB ZRPE3783 A1 = RB-RA*(D/C) ZRPE3784 A7 = RA+G*D+H*F ZRPE3785 RETURN ZRPE3786 END ZRPE3787 SUBROUTINE ZRPQLF (ITYPE) ZRPF3816 INTEGER ITYPE ZRPF3818 INTEGER N,NN,I ZRPF3820 REAL*8 ARE,ETA,RMRE ZRPF3821 DOUBLE PRECISION P(101),QP(101),RK(101),QK(101),SVK(101) ZRPF3822 DOUBLE PRECISION SR,SI,U,V,RA,RB,C,D,A1,A2,A3, ZRPF3823 1 A6,A7,E,F,G,H,SZR,SZI,RLZR,RLZI,TEMP,ZERO ZRPF3824 COMMON /ZRPQLJ/ P,QP,RK,QK,SVK,SR,SI,U,V,RA,RB,C,D,A1,A2,A3,A6,ZRPF3825 1 A7,E,F,G,H,SZR,SZI,RLZR,RLZI,ETA,ARE,RMRE,N,NN ZRPF3826 DATA ZERO/0.0D0/ ZRPF3827 IF (ITYPE.EQ.3) GO TO 20 ZRPF3831 TEMP = RA ZRPF3832 IF (ITYPE.EQ.1) TEMP = RB ZRPF3833 IF (DABS(A1).GT.DABS(TEMP)*ETA*10.D0) GO TO 10 ZRPF3834 RK(1) = ZERO ZRPF3837 RK(2) = -A7*QP(1) ZRPF3838 DO 5 I=3,N ZRPF3839 5 RK(I) = A3*QK(I-2)-A7*QP(I-1) ZRPF3840 RETURN ZRPF3841 10 A7 = A7/A1 ZRPF3843 A3 = A3/A1 ZRPF3844 RK(1) = QP(1) ZRPF3845 RK(2) = QP(2)-A7*QP(1) ZRPF3846 DO 15 I=3,N ZRPF3847 15 RK(I) = A3*QK(I-2)-A7*QP(I-1)+QP(I) ZRPF3848 RETURN ZRPF3849 20 RK(1) = ZERO ZRPF3852 RK(2) = ZERO ZRPF3853 DO 25 I=3,N ZRPF3854 25 RK(I) = QK(I-2) ZRPF3855 RETURN ZRPF3856 END ZRPF3857 SUBROUTINE ZRPQLG (ITYPE,UU,VV) ZRPG3886 INTEGER ITYPE ZRPG3888 DOUBLE PRECISION UU,VV ZRPG3889 INTEGER N,NN ZRPG3891 REAL*8 ARE,ETA,RMRE ZRPG3892 DOUBLE PRECISION P(101),QP(101),RK(101),QK(101),SVK(101) ZRPG3893 DOUBLE PRECISION SR,SI,U,V,RA,RB,C,D,A1,A2,A3, ZRPG3894 1 A6,A7,E,F,G,H,SZR,SZI,RLZR,RLZI, ZRPG3895 2 A4,A5,B1,B2,C1,C2,C3,C4,TEMP,ZERO ZRPG3896 COMMON /ZRPQLJ/ P,QP,RK,QK,SVK,SR,SI,U,V,RA,RB,C,D,A1,A2,A3,A6,ZRPG3897 1 A7,E,F,G,H,SZR,SZI,RLZR,RLZI,ETA,ARE,RMRE,N,NN ZRPG3898 DATA ZERO/0.0D0/ ZRPG3899 IF (ITYPE.EQ.3) GO TO 15 ZRPG3906 IF (ITYPE.EQ.2) GO TO 5 ZRPG3907 A4 = RA+U*RB+H*F ZRPG3908 A5 = C+(U+V*F)*D ZRPG3909 GO TO 10 ZRPG3910 5 A4 = (RA+G)*F+H ZRPG3911 A5 = (F+U)*C+V*D ZRPG3912 10 B1 = -RK(N)/P(NN) ZRPG3915 B2 = -(RK(N-1)+B1*P(N))/P(NN) ZRPG3916 C1 = V*B2*A1 ZRPG3917 C2 = B1*A7 ZRPG3918 C3 = B1*B1*A3 ZRPG3919 C4 = C1-C2-C3 ZRPG3920 TEMP = A5+B1*A4-C4 ZRPG3921 IF (TEMP.EQ.ZERO) GO TO 15 ZRPG3922 UU = U-(U*(C3+C2)+V*(B1*A1+B2*A7))/TEMP ZRPG3923 VV = V*(1+C4/TEMP) ZRPG3924 RETURN ZRPG3925 15 UU = ZERO ZRPG3927 VV = ZERO ZRPG3928 RETURN ZRPG3929 END ZRPG3930 SUBROUTINE ZRPQLH (NN,U,V,P,Q,RA,RB) ZRPH3959 INTEGER NN ZRPH3961 DOUBLE PRECISION P(NN),Q(NN),U,V,RA,RB ZRPH3962 INTEGER I ZRPH3964 DOUBLE PRECISION C ZRPH3965 RB = P(1) ZRPH3970 Q(1) = RB ZRPH3971 RA = P(2)-U*RB ZRPH3972 Q(2) = RA ZRPH3973 DO 5 I=3,NN ZRPH3974 C = P(I)-U*RA-V*RB ZRPH3975 Q(I) = C ZRPH3976 RB = RA ZRPH3977 RA = C ZRPH3978 5 CONTINUE ZRPH3979 RETURN ZRPH3980 END ZRPH3981 SUBROUTINE ZRPQLI (RA,B1,C,SR,SI,RLR,RLI) ZRPI4010 DOUBLE PRECISION RA,B1,C,SR,SI,RLR,RLI ZRPI4012 DOUBLE PRECISION RB,D,E,ZERO,ONE,TWO ZRPI4014 DATA ZERO,ONE,TWO/0.0D0,1.0D0,2.0D0/ ZRPI4015 IF (RA.NE.ZERO) GO TO 10 ZRPI4026 SR = ZERO ZRPI4027 IF (B1.NE.ZERO) SR = -C/B1 ZRPI4028 RLR = ZERO ZRPI4029 5 SI = ZERO ZRPI4030 RLI = ZERO ZRPI4031 RETURN ZRPI4032 10 IF (C.NE.ZERO) GO TO 15 ZRPI4033 SR = ZERO ZRPI4034 RLR = -B1/RA ZRPI4035 GO TO 5 ZRPI4036 15 RB = B1/TWO ZRPI4039 IF (DABS(RB).LT.DABS(C)) GO TO 20 ZR 4040 E = ONE-(RA/RB)*(C/RB) ZRPI4041 D = DSQRT(DABS(E))*DABS(RB) ZR 4042 GO TO 25 ZRPI4043 20 E = RA ZRPI4044 IF (C.LT.ZERO) E = -RA ZRPI4045 E = RB*(RB/DABS(C))-E ZRP 4046 D = DSQRT(DABS(E))*DSQRT(DABS(C)) ZR 4047 25 IF (E.LT.ZERO) GO TO 30 ZRPI4048 IF (RB.GE.ZERO) D = -D ZRPI4050 RLR = (-RB+D)/RA ZRPI4051 SR = ZERO ZRPI4052 IF (RLR.NE.ZERO) SR = (C/RLR)/RA ZRPI4053 GO TO 5 ZRPI4054 30 SR = -RB/RA ZRPI4056 RLR = SR ZRPI4057 SI = DABS(D/RA) ZRP 4058 RLI = -SI ZRPI4059 RETURN ZRPI4060 END ZRPI4061 C********************************************************************** C C SUBROUTINE PPLOT C C PROVIDES ROOT LOCUS PLOTS C********************************************************************** SUBROUTINE PPLOT(KP,NR) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*1 ICH,IPR COMMON/LOCO/XX(200,10),YY(200,10) DIMENSION NX(200,9),NY(200,9),IPR(106),ICH(12) DATA ICH/'1','2','3','4','5','6','7','8','9',' ','-','|'/ IER=0 IF(KP.GT.200)KP=200 IF(NR.GT.9)NR=9 JW=3 ZERO=0.0D0 AX=35.0D0 AY=28.0D0 XMIN=0.0D0 YMIN=0.0D0 XMAX=0.0D0 YMAX=0.0D0 C C DETERMINE SCALE FACTORS C DO 20 J=1,NR DO 20 I=1,KP IF(XX(I,J).LT.XMIN)XMIN=XX(I,J) IF(XX(I,J).GT.XMAX)XMAX=XX(I,J) IF(YY(I,J).LT.YMIN)YMIN=YY(I,J) IF(YY(I,J).GT.YMAX)YMAX=YY(I,J) 20 CONTINUE IF(XMIN.LT.-71.0D0) XMIN=-71.0D0 IF(XMIN.EQ.0.0D0)XMIN=-1.0D0 IF(XMAX.EQ.0.0D0)XMAX=-XMIN/2.0D0 IF(YMAX.EQ.0.0D0)YMAX=XMAX IF(YMIN.EQ.0.0D0)YMIN=-YMAX RX=-XMIN/XMAX IF(RX.GT.2.0D0)XMAX=-XMIN/2.0D0 RY=-YMAX/YMIN IF(RY.LT.1.0D0)YMAX=-YMIN RY=XMAX/YMAX IF(RY.LT.1.0D0)XMAX=YMAX IF(RY.GT.1.0D0)YMAX=XMAX YMIN=-YMAX XMIN=-2.0D0*XMAX AX=AX/XMAX AY=AY/YMAX C C SCALE DATA AND PROTECT AGAINST ARRAY OVERFLOW C DO 30 J=1,NR DO 30 I=1,KP XX(I,J)=AX*(XX(I,J)-XMIN)+0.1D0 YY(I,J)=AY*(YY(I,J)-YMIN)+0.1D0 NX(I,J)=IDINT(XX(I,J))+1 IF(NX(I,J).LT.2)NX(I,J)=2 IF(NX(I,J).GT.105) NX(I,J)=105 NY(I,J)=57-IDINT(YY(I,J)) IF(NY(I,J).LT.2)NY(I,J)=2 IF(NY(I,J).GT.56)NY(I,J)=56 30 CONTINUE C C CONSTRUCT PLOT C IPR(1)=ICH(10) IPR(106)=ICH(10) KJ=1 DO 40 K=2,105 40 IPR(K)=ICH(11) WRITE(JW,200)IPR,KJ 200 FORMAT(10X,106A1,12X,I2) IPR(1)=ICH(12) IPR(106)=ICH(12) KJ=2 100 IF(KJ.NE.29)GO TO 60 DO 50 K=2,105 50 IPR(K)=ICH(11) IPR(101)=ICH(11) GO TO 80 60 DO 70 K=2,105 70 IPR(K)=ICH(10) IPR(71)=ICH(12) C C LOAD IPR ARRAY WITH APPROPRIATE NUMBERS C 80 DO 90 J=1,NR DO 90 K=1,KP 90 IF(NY(K,J).EQ.KJ)IPR(NX(K,J))=ICH(J) IPR(1)=ICH(12) IPR(106)=ICH(12) WRITE(JW,200)IPR,KJ IF(KJ.EQ.29)WRITE(JW,210)ZERO KJ=KJ+1 IF(KJ.GT.56) GO TO 101 GO TO 100 210 FORMAT(1H+,3X,F3.1) 101 IPR(1)=ICH(10) IPR(106)=ICH(10) DO 110 K=2,105 110 IPR(K)=ICH(11) WRITE(JW,200)IPR,KJ WRITE(JW,220)ZERO 220 FORMAT(79X,F3.1) C C WRITE OUT MAXIMUM AND MINIMUM VALUES PLOTTED C WRITE(JW,230)XMIN,XMAX 230 FORMAT(13X,'XMIN=',F12.4,3X,'XMAX=',F12.4,3X) WRITE(JW,231) YMIN,YMAX 231 FORMAT(13X,'YMIN=',F12.4,3X,'YMAX=',F12.4) RETURN END C********************************************************************** C C SUBROUTINE UERTST C C PROVIDES ERROR MESSAGES C********************************************************************** SUBROUTINE UERTST(IER,NAME) CHARACTER*4 ITYP DIMENSION ITYP(5,4),IBIT(4) INTEGER*2 NAME(3) INTEGER WARN,WARF,TERM,JWRITE EQUIVALENCE (IBIT(1),WARN),(IBIT(2),WARF),(IBIT(3),TERM) DATA ITYP/'WARN','ING ',' ',' ',' ','WARN','ING(','WITH', 1' FIX',') ','TERM','INAL',' ',' ',' ','NON-','DEFI','NE 1D ',' ',' '/,IBIT/32,64,128,0/ C C IDENTIFY WRITE UNIT C JWRITE=3 IER2=IER IF (IER2.GE.WARN) GO TO 1 C*** UNDEFINED IER1=4 GO TO 3 1 IF (IER2.LT.TERM) GO TO 2 C*** TERMINAL IER1=3 GO TO 3 2 IF (IER2.LT.WARF) GO TO 4 C*** WARNING (WITH FIX) IER1=2 GO TO 3 C*** WARNING C*** EXTRACT 'N' 3 IER2=IER2-IBIT(IER1) 4 IER1=1 C*** PRINT ERROR MESSAGE WRITE (JWRITE,5) (ITYP(I,IER1),I=1,5),NAME,IER2,IER 5 FORMAT (1X,/,1X,'ERROR MESSAGE FROM UERTST',2X,5A4,4X,3A2,4X,I2,2X 1,'IER= ',I3) RETURN END C************************************************************************ 4955 C C 4956 C C 4957 C-----------------------------------------------------------------------C 4958 C SUBROUTINE GENPT(NCSU VERSION) C 4959 C-----------------------------------------------------------------------C 4960 C C 4961 C C 4962 C THIS IS THE STANDARD PRINTER PLOT ROUTINE USED C 4963 C AT NCSU. C 4964 C C 4965 C C 4966 C C 4967 C THE SOURCE OF THIS ROUTINE IS UNKNOWN TO PRESENT C 4968 C COMPUTER CENTER STAFF MEMBERS, C 4969 C C 4970 C C 4971 C NOTE: ALTHOUGH 200 POINTS MAY BE PASSED TO THIS C 4972 C VERSION OF GENPT, ONLY 100 POINTS CAN BE C 4973 C PLOTTED SINCE THE GRAPH IS 100 SPACES BY 50 C 4974 C LINES. THIS ROUTINE DOES NOT PROVIDE THE C 4975 C CAPABILITY TO LABEL THE AXES. C 4976 C C 4977 C C 4978 C VERSION COPIED: JUNE 1981 C 4979 C FORTRAN 77 VERSION PREPARED OCTOBER 1981 C 4980 C C 4981 C C 4982 C************************************************************************ 4983 C C SUBROUTINE GENPT(X,Y,NOPTS,ITERM) REAL X(NOPTS),Y(NOPTS) REAL TEMP,DT(2),XP(11) CHARACTER*1 CA,BLK,A(101) INTEGER*4 NNCT DATA BLK,CA/' ','.'/ C IF(ITERM.EQ.0)GO TO 1 DO 200 I=1,NOPTS TEMP=X(I) X(I)=Y(I) Y(I)=TEMP 200 CONTINUE 1 XMIN=X(1) XMAX=X(1) YMIN=Y(1) YMAX=Y(1) DO 3 I=2,NOPTS XMAX=AMAX1(XMAX,X(I)) YMAX=AMAX1(YMAX,Y(I)) XMIN=AMIN1(XMIN,X(I)) 3 YMIN=AMIN1(YMIN,Y(I)) NNCT=100 C CALL SCALE2 (XMIN,XMAX,NNCT,BOTX,TOPX,NAL) C NINT=100 IF(ITERM.NE.0) NINT=50 DT(1)=(TOPX-BOTX)/NINT C CALL SCALE2 (YMIN,YMAX,NNCT,BOTY,TOPY,NAL) C NINT=50 IF(ITERM.EQ.1) NINT=100 DT(2)=(TOPY-BOTY)/NINT IF(ITERM.EQ.0)GO TO 202 IL=101 IF(ITERM.EQ.2) IL=51 IM=IL/10 DXSC=5.*DT(1) DO 217 I=1,11 217 XP(I)=BOTX+FLOAT(I-1)*DXSC WRITE (3,199) WRITE (3,111) DY2=0.500001*DT(2) DO 222 II=1,IL I=II-1 DO 265 J=1,51 265 A(J)=BLK YN=TOPY-FLOAT(I)*DT(2) DO 218 J=1,NOPTS IF(ABS(YN-Y(J)).GT.DY2)GO TO 218 L=1+IFIX(.5+(X(J)-BOTX)/DT(1)) A(L)=CA 218 CONTINUE IF(MOD(I,IM).EQ.0)GO TO 220 WRITE (3,112) (A(J),J=1,51) GO TO 222 220 WRITE (3,113) YN,(A(J),J=1,51) 222 CONTINUE WRITE (3,114)(XP(I),I=1,11) WRITE (3,115) XMAX,XMIN,YMAX,YMIN,DT(1),DT(2) DO 250 I=1,NOPTS TEMP=X(I) X(I)=Y(I) Y(I)=TEMP 250 CONTINUE RETURN 202 CONTINUE DXSC=10.*DT(1) DO 17 I=1,11 17 XP(I)=BOTX+FLOAT(I-1)*DXSC WRITE(3,101) DY2=0.500001*DT(2) DO 22 II=1,51 I=II-1 DO 165 J=1,101 165 A(J)=BLK YN=TOPY-FLOAT(I)*DT(2) DO 18 J=1,NOPTS IF(ABS(YN-Y(J)).GT.DY2) GO TO 18 L=1+IFIX(.5+(X(J)-BOTX)/DT(1)) A(L)=CA 18 CONTINUE IF(MOD(I,5).EQ.0) GO TO 20 WRITE (3,102) (A(J),J=1,101) GO TO 22 20 WRITE (3,103) YN,(A(J),J=1,101) 22 CONTINUE WRITE (3,104) (XP(I),I=1,11) WRITE (3,105) XMAX,XMIN,YMAX,YMIN,DT(1),DT(2) RETURN 101 FORMAT(1H1,14X,1H*,10(9X,1H*)/14X,103(1H*)) 199 FORMAT(1H1,29X,21HX AND Y AXES SWITCHED ) 111 FORMAT(15X,1H*,10(4X,1H*)/14X,53(1H*)) 102 FORMAT(14X,1H*,101A1,1H*) 112 FORMAT(14X,1H*,51A1,1H*) 103 FORMAT(1X,1PE11.2,3H **,101A1,2H**) 113 FORMAT(1X,1PE11.2,3H **,51A1,2H**) 104 FORMAT(14X,103(1H*)/6X,11(9X,1H*)/9X,11(1PE10.2)/) 114 FORMAT(14X,53(1H*)/11X,11(4X,1H*)/11X,6HY-AXIS,11(/9X,1PE10.2)/) 105 FORMAT(1X,5HXMAX=,1PE14.5,2X,5HXMIN=,E14.5/2X,5HYMAX=,E14.5,2X, 15HYMIN=,E14.5/2X,6HXINCR=,E14.5,2X,6HYINCR=,E14.5) 115 FORMAT(1X,5HYMAX=,1PE14.5,2X,5HYMIN=,E14.5/2X,5HXMAX=,E14.5,2X, 15HXMIN=,E14.5/2X,6HYINCR=,E14.5,2X,6HXINCR=,E14.5) END C********************************************************************** C C SUBROUTINE SCAL2 C C********************************************************************** SUBROUTINE SCALE2 (DMIN,DMAX,N,SMIN,SMAX,NAL) INTEGER*4 N DIMENSION VINT(20) DATA VINT/1.,2.,3.,4.,5.,6.,7.,8.,9.,0.,0.,0.,0.,0.,0.,0.,0.,0. 1,0.,0./ DATA NORND/9/ DATA DEL/2.E-6/ C A=(DMAX-DMIN)/FLOAT(N) NAL=ALOG10(A) IF (A.LT.1.) NAL=NAL-1 B=A/(10.**NAL) DO 20 I=1,NORND IF (B.LT.(VINT(I)+DEL)) GO TO 30 20 CONTINUE 25 I=1 NAL=NAL+1 30 DIST=VINT(I)*10.**NAL FM=DMIN/DIST M1=FM IF(FM.LT.0.) M1=M1-1 IF(IABS(M1)+1.-FM.LT.DEL) M1=M1+1 FM=DMAX/DIST M2=FM+1. IF (FM.LT.-1) M2=M2-1 IF (ABS(FM+1.-M2).LT.DEL) M2=M2-1 NP=M2-M1 IF (NP.LE.N) GO TO 40 I=I+1 IF (I-NORND)30,30,25 40 NP=(N-NP)/2 SMIN=(M1-NP)*DIST 7 SMAX=SMIN+N*DIST IF (SMIN.GT.DMIN) SMIN=DMIN IF (SMAX.LT.DMAX) SMAX=DMAX RETURN END C***********************************************************************C 7985 C C 7986 C-----------------------------------------------------------------------C 7987 C SUBROUTINE NEWTON C 7988 C-----------------------------------------------------------------------C 7989 C C 7990 C PURPOSE: TO REFINE ROOT VALUES OBTAINED BY ZRPOLY. C 7991 C METHOD: NEWTON-RAPHSON IN EXTENDED PRECISION C 7992 C VARIABLE IDENTIFICATION: C 7993 C C 7994 C BZ - COEFFICIENTS OF CHARACTERISTIC EQUATION C 7995 C ZF - REAL AND IMAGINARY PARTS OD ROOTS OF C.E. C 7996 C BZP - COEFFICIENTS OF DERIVATIVE OF C.E. C 7997 C RES - RESIDUAL OF C.E. WITH ROOT VALUES SUBSTITUTED C 7998 C OR RESIDUAL OF DERIVATIVE OF C.E. C 7999 C MX - NUMBER OF COEFFICIENTS IN C.E. C 8000 C MQ - STORAGE LOCATION OF SPECIFIC IMAGINARY C 8001 C PART OF ROOT C 8002 C C 8003 C DATE OF CURRENT VERSION: JUNE 1985 C 8004 C C 8005 C***********************************************************************C 8006 C C 8007 SUBROUTINE NEWTON(BZ,ZF,MQ,MX) 8008 IMPLICIT REAL*8(A-H,O-Z) 8009 REAL*8 RES,RNUM 8010 DIMENSION BZ(11),ZF(20),BZP(10) 8011 C 8012 MR=MQ-2 8013 MX1=MX-1 8014 DO 1 J=1,MX1 8015 JJ=MX1-J+1 8016 1 BZP(J)=BZ(J)*DBLE(FLOAT(JJ)) 8017 DO 4 I=1,2 8018 II=I*2 8019 LCNT=0 8020 5 RES=0.0D0 8021 DO 2 J=1,MX 8022 2 RES=RES*ZF(MR+II)+BZ(J) 8023 IF(DABS(RES).LT.1.0D-50) GO TO 4 8024 IF(LCNT.GT.100) GO TO 4 8025 RNUM=RES 8026 RES=0.0D0 8027 DO 3 J=1,MX1 8028 3 RES=RES*ZF(MR+II)+BZP(J) 8029 ZF(MR+II)=ZF(MR+II)-RNUM/RES 8030 LCNT=LCNT+1 8031 GO TO 5 8032 4 CONTINUE 8033 DO 7 J=1,8,4 8034 DOFF=ZF(J)-ZF(J+2) 8035 7 IF(DABS(DOFF).LT.1.0D-15) ZF(J)=ZF(J+2) 8036 RETURN 8037 END 8038