1C $Header: /var/cvs/mbdyn/mbdyn/mbdyn-1.0/mbdyn/aero/aerod2.f,v 1.31 2017/01/12 14:45:58 masarati Exp $ 2C MBDyn (C) is a multibody analysis code. 3C http://www.mbdyn.org 4C 5C Copyright (C) 1996-2017 6C 7C Pierangelo Masarati <masarati@aero.polimi.it> 8C Paolo Mantegazza <mantegazza@aero.polimi.it> 9C 10C Dipartimento di Ingegneria Aerospaziale - Politecnico di Milano 11C via La Masa, 34 - 20156 Milano, Italy 12C http://www.aero.polimi.it 13C 14C Changing this copyright notice is forbidden. 15C 16C This program is free software; you can redistribute it and/or modify 17C it under the terms of the GNU General Public License as published by 18C the Free Software Foundation (version 2 of the License). 19C 20C 21C This program is distributed in the hope that it will be useful, 22C but WITHOUT ANY WARRANTY; without even the implied warranty of 23C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 24C GNU General Public License for more details. 25C 26C You should have received a copy of the GNU General Public License 27C along with this program; if not, write to the Free Software 28C Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 29C 30C ============================================================================ 31C 32C These routines have been supplied by 33C Professor Massimiliano Lanz <massimiliano.lanz@polimi.it> 34C 35C ============================================================================ 36C 37C************************* AEROD ********************* 38c MODIFICATA NEL CALCOLO TNG 39C= COMPILER (LINK=IBJ$) 40 SUBROUTINE AEROD2(W, VAM, TNG, OUTA, INST, RSPEED, JPRO) 41C 42C W: vettore lungo 6; velocita' del corpo aerodinamico nel sistema locale, 43C riferita ad un punto arbitrario 44C VAM: vettore lungo 6, contiene dati del problema 45C TNG: vettore lungo 12, presumubilmente il termine noto 46C (ora l'ho ridotto a 6) 47C OUTA: vettore lungo 20, usato in i/o con subroutines di aerod 48C INST: flag di instazionario (>0) 49C RSPEED -> MODOMEGA: modulo della velocita' di rotazione 50C JPRO: tipo di profilo 51 52 53 IMPLICIT NONE 54 55C Input/Output 56 REAL*8 W(6), VAM(6), TNG(6), OUTA(*), RSPEED 57 INTEGER*4 INST, JPRO 58 59C Local 60 REAL*8 DENS, CS, CORDA, B, D, VP2, VP, V, RK 61 INTEGER*4 IGO, I 62 63 REAL*8 CLIFT, CDRAG, CMOME, ASLOP, BSLOP, 64 & CSLOP, CRF, CFSLOP, DCPDM, DCDRDM, DCMDM, DCRFDM 65 66 REAL*8 DM(3,6), VCSTR(6), BM(6,4), AASTR(6,6), 67 & A1STR(4,3), AP(4,6), QASTR(6) 68C 69C DEFINIZIONI VETTORE VAM 70C 71C VAM(1): densita' dell'aria 72C VAM(2): celerita' del suono 73C VAM(3): corda 74C VAM(4): 1/4 corda 75C VAM(5): 3/4 corda 76C VAM(6): svergolamento (non usato; aggiunto esternamente alla rotaz. profilo) 77C 78 DENS = VAM(1) 79 CS = VAM(2) 80 CORDA = VAM(3) 81 B = VAM(4) 82 D = VAM(5) 83C SVER = VAM(6) 84C 85C COSTRUZIONE VC* 86C 87 CALL DZERO(DM, 18) 88 DO 10 I = 1,3 89 DM(I,I) = 1.D0 90 10 CONTINUE 91 DM(2,6) = D 92CCC Ci vuole anche questo ?!? 93 DM(3,5) = -D 94C 95C Ora DM ha la struttura seguente: 96C [ 1 0 0 0 0 0 ] 97C [ 0 1 0 0 0 D ] 98C [ 0 0 1 0 -D 0 ] 99C 100C e quindi se la velocita' W del corpo aerodinamico e' organizzata come 101C W = [v1 v2 v3 w1 w2 w3 ] 102C con: 103C 1 direzione di avanzamento, 104C 2 direzione normale, 105C 3 direzione lungo l'apertura, 106C DM*W da' la velocita' nel punto a 3/4 della corda 107C 108C Moltiplica DM per W a dare VCSTR (velocita' nel sistema locale, punto 3/4 c) 109 CALL DPROMV(DM, 3, 3, 6, W, VCSTR, 0) 110C 111C VP2 e' il modulo della velocita' nel piano del profilo al quadrato 112C VP e' il modulo della velocita' nel piano del profilo 113C V e' il modulo della velocita' 114 VP2 = VCSTR(1)**2+VCSTR(2)**2 115 VP = DSQRT(VP2) 116 V = DSQRT(VP2+VCSTR(3)**2) 117C 118C Se il numero di Mach e' troppo piccolo, si ferma 119 IF(V/CS.LT.1.D-6) THEN 120 CALL DZERO(TNG, 6) 121 RETURN 122 ENDIF 123C 124C VP: modulo della velocita' nel piano del profilo 125C V: modulo della velocita' totale 126C 127C CALCOLO COEFFICIENTI AERODINAMICI 128C 129C INST e' il flag di instazionarieta'; 130C puo' valere 0, 1, 2 e forza l'utilizzo del metodo PK del relativo ordine 131C FIXME: solo 0 e' valido, gli altri due sono probabilmente bacati ... 132 IGO = INST+1 133 GOTO(100, 200, 300), IGO 134C 135C Steady 136 100 CALL COE0(VCSTR, OUTA, CS, 137 & CLIFT, CDRAG, CMOME, ASLOP, BSLOP, CSLOP, 138 & CRF, CFSLOP, DCPDM, DCDRDM, DCMDM, DCRFDM, JPRO) 139 GOTO 400 140C 141C Unsteady, HARRIS A.H.S. JULY 1970 142 200 CALL COE1(VCSTR, OUTA, CS, CORDA, RSPEED, 143 & CLIFT, CDRAG, CMOME, ASLOP, BSLOP, CSLOP, 144 & CRF, CFSLOP, DCPDM, DCDRDM, DCMDM, DCRFDM, JPRO) 145 GOTO 400 146C 147C Unsteady, BIELAWA 31TH A.H.S. FORUM 1975 148 300 CALL COE2(VCSTR, OUTA, CS, CORDA, RSPEED, 149 & CLIFT, CDRAG, CMOME, ASLOP, BSLOP, CSLOP, 150 & CRF, CFSLOP, DCPDM, DCDRDM, DCMDM, DCRFDM, JPRO) 151 152C 153C COSTRUZIONE A1* ? 154C 155 400 CONTINUE 156 CALL DZERO(A1STR, 12) 157 RK = -.5D0*DENS*CORDA 158 A1STR(1,1) = RK*(CRF*(VP-V)-CDRAG*VP) 159 A1STR(1,2) = -RK*CLIFT*VP 160 A1STR(2,1) = -A1STR(1, 2) 161 A1STR(2,2) = A1STR(1, 1) 162 A1STR(3,3) = -RK*CRF*V 163 A1STR(4,1) = -RK*CMOME*CORDA*VCSTR(1) 164 A1STR(4,2) = -RK*CMOME*CORDA*VCSTR(2) 165C 166C La matrice BM ha la forma 167C 168C [ 1 0 0 0 ] 169C [ 0 1 0 0 ] 170C [ 0 0 1 0 ] 171C 172C [ 0 0 0 0 ] 173C [ 0 0 0 0 ] 174C [ 0 B 0 1 ] 175C 176C dove B e' la posizione del punto a 1/4 della corda rispetto al riferimento 177C (punto di applicazione delle forze aerodinamiche) 178C 179 CALL DZERO(BM, 24) 180 DO 20 I = 1,3 181 BM(I,I) = 1.D0 182 20 CONTINUE 183 BM(6,2) = B 184 BM(6,4) = 1.D0 185 CALL DPROMM(A1STR, 4, 4, 3, DM, 3, AP, 4, 6, 0, 0) 186 CALL DPROMM(BM, 6, 6, 4, AP, 4, AASTR, 6, 6, 0, 0) 187C 188C COSTRUZIONE Q* ? 189C 190CCCCC Copia gli ultimi 3 termini di W in VCSTR (elimina l'effetto della 191CCCCC moltiplicazione iniziale per DM) 192 CALL MOVE(VCSTR(4), W(4), 3) 193CCCCC Quindi moltiplica la matrice AASTR per VCSTR a dare QASTR, le forze 194 CALL DPROMV(AASTR, 6, 6, 6, VCSTR, QASTR, 0) 195C cambio segno alle forze perch� le calcola con segno opposto 196 do I = 1,6 197 TNG(I) = -QASTR(I) 198 end do 199C 200 END 201 202 203 204C************************* COE0 ********************** 205C= COMPILER (LINK=IBJ$) 206 SUBROUTINE COE0(VCSTR, OUTA, CS, 207 & CLIFT, CDRAG, CMOME, ASLOP, BSLOP, CSLOP, 208 & CRF, CFSLOP, DCPDM, DCRDRM, DCMDM, DCRFDM, JPRO) 209 210 211C Modificato da Pierangelo Masarati 19/11/97 212 implicit none 213C 214C Input 215 real*8 VCSTR(*), CS 216 integer*4 JPRO 217C 218C Output 219 real*8 OUTA(*), 220 & CLIFT,CDRAG,CMOME,ASLOP,BSLOP,CSLOP,CRF,CFSLOP,DCPDM,DCRDRM, 221 & DCMDM,DCRFDM 222C 223C Local 224 real*8 ALFA,VC1,GAM,COSGAM,RMACH,ASLRF,ASLOP0,CSLOP0 225C 226C Data 227 real*8 PG, DEGRAD 228 DATA PG /3.14159265358979323846D0/, DEGRAD /.017453293D0/ 229C 230C CALCOLA I COEFFICIENTI AERODINAMICI CON TEORIA STAZIONARIA 231C (CON CORREZIONE PER FRECCIA SECONDO HARRIS A.H.S.JULY 1970) 232C 233 ALFA = DATAN2(-VCSTR(2),VCSTR(1)) 234 OUTA(2) = ALFA/DEGRAD 235 VC1 = DABS(VCSTR(1)) 236 GAM = DATAN2(-VCSTR(3),VC1) 237 OUTA(3) = GAM/DEGRAD 238 IF(DABS(GAM).GT.PG/3.D0) GAM = PG/3.D0 239 COSGAM = DCOS(GAM) 240 RMACH = DSQRT(VCSTR(1)**2+VCSTR(2)**2+VCSTR(3)**2)/CS 241 RMACH = RMACH*DSQRT(COSGAM) 242 OUTA(4) = RMACH 243 CALL CPCRCM(ALFA, RMACH, 244 & CLIFT, CDRAG, CMOME, ASLOP, BSLOP, CSLOP, 245 & CRF, CFSLOP, DCPDM, DCRDRM, DCMDM, DCRFDM, 246 & ASLOP0, CSLOP0, JPRO, 0) 247 ASLRF = ASLOP0 248 IF(DABS(ALFA).LT.1.D-6) GOTO 10 249 ASLRF = CLIFT/(ALFA*COSGAM) 250 IF(ASLRF.GT.ASLOP0) ASLRF = ASLOP0 251 10 CLIFT = ASLRF*ALFA 252 OUTA(5) = CLIFT 253 OUTA(6) = CDRAG 254 OUTA(7) = CMOME 255C 256 RETURN 257C 258C 259 END 260C************************* COEPRD ********************** 261C= COMPILER (LINK=IBJ$) 262 SUBROUTINE COEPRD(DA, OUTA) 263 IMPLICIT NONE 264 REAL*8 DA, OUTA(*) 265 REAL*8 ALF1, ALF2 266 ALF1 = OUTA(9)/DA 267 OUTA(9) = ALF1 268 ALF2 = OUTA(10)/(DA*DA) 269 OUTA(10) = ALF2 270 RETURN 271 END 272C************************* COE1 ********************** 273C= COMPILER (LINK=IBJ$) 274 SUBROUTINE COE1(VCSTR, OUTA, CS, CORDA, RSPEED, 275 & CLIFT, CDRAG, CMOME, ASLOP, BSLOP, CSLOP, 276 & CRF, CFSLOP, DCPDM, DCRDRM, DCMDM, DCRFDM, JPRO) 277 278 279 implicit none 280C IMPLICIT REAL*8(A-H,O-Z) 281C 282C Input: 283 real*8 VCSTR(*),CS,CORDA,RSPEED 284 integer*4 JPRO 285C 286C Output: 287 real*8 OUTA(*), 288 & CLIFT,CDRAG,CMOME,ASLOP,BSLOP,CSLOP,CRF,CFSLOP, 289 & DCPDM,DCRDRM,DCMDM,DCRFDM 290C 291C Functions: 292 real*8 THF,THG 293C 294C Local: 295 real*8 ABE,VC1,ALF1,ALF2,VP,GAM,COSGAM,RMACH,FR,RMM,DALF, 296 & SEGNO,AREF,ASLOP0,CSLOP0,ASLRF,RK,PA,ALF 297C 298C Data: 299 real*8 PG,DEGRAD 300 DATA PG /3.14159265358979323846D0/, DEGRAD /.017453293D0/ 301C 302C CALCOLA COEFFICIENTI AERODINAMICI CON TEORIA INSTAZIONARIA 303C SECONDO TEORIA HARRIS A.H.S. JULY 1970 304C 305 ABE = DATAN2(-VCSTR(2),VCSTR(1)) 306 OUTA(2) = ABE/DEGRAD 307 VC1 = DABS(VCSTR(1)) 308C 309C Questa operazione e' stata spostata all'esterno (COEPRD) per consentire 310C l'uso di integratori impliciti, per i quali la correzione deve essere fatta 311C solo a processo iterativo concluso 312C 313C DAA = DA*RSPEED 314C ALF1 = OUTA(9)/DAA 315C OUTA(9) = ALF1 316C ALF2 = OUTA(10)/(DAA*DAA) 317C OUTA(10) = ALF2 318C 319C Sono: ALF1 = theta/t 320C ALF2 = theta/tt 321 ALF1 = OUTA(9) 322 ALF2 = OUTA(10) 323C 324 VP = DSQRT(VCSTR(1)**2+VCSTR(2)**2) 325 GAM = DATAN2(-VCSTR(3),VC1) 326 OUTA(3) = GAM/DEGRAD 327 IF(DABS(GAM).GT.PG/3.D0) GAM = PG/3.D0 328 COSGAM = DCOS(GAM) 329 RMACH = DSQRT(VCSTR(1)**2+VCSTR(2)**2+VCSTR(3)**2)/CS 330 RMACH = RMACH*DSQRT(COSGAM) 331 OUTA(4) = RMACH 332C Viene usato ALF1 per la prima volta: correzione instazionaria 333 FR = .5D0*CORDA*ALF1*RSPEED/VP 334 FR = DABS(FR) 335 OUTA(11) = FR 336C Mach effettivo usato solo tra .3 e .6 337 RMM = .3D0 338 IF(RMACH.GT..3D0) RMM = RMACH 339 IF(RMACH.GT..6D0) RMM = .6D0 340C 341 DALF = 61.5D0*DLOG(.6D0/RMM)*DSQRT(FR)*PG/180.D0 342 SEGNO = 1.D0 343 IF(ALF1.LT.0.D0) SEGNO = -1.D0 344 AREF = ABE-SEGNO*DALF 345 OUTA(12) = AREF/DEGRAD 346 CALL CPCRCM(AREF, RMACH, 347 & CLIFT, CDRAG, CMOME, ASLOP, BSLOP, CSLOP, 348 & CRF, CFSLOP, DCPDM, DCRDRM, DCMDM, DCRFDM, 349 & ASLOP0, CSLOP0, JPRO, 0) 350 ASLRF = ASLOP0 351C Se l'angolo di incidenza e' significativo corregge la slope di CL/alfa 352 IF(DABS(AREF).LT.1.D-6) GOTO 10 353 ASLRF = CLIFT/(AREF*COSGAM) 354 IF(ASLRF.GT.ASLOP0) ASLRF = ASLOP0 355 10 CONTINUE 356C Frequenza ridotta (spero che con VP = 0 non si entri qui ...) 357 RK = .5D0*CORDA*RSPEED/VP 358 PA = .25D0 359C Corregge l'angolo di incidenza attraverso le funzioni di Theodorsen 360 ALF = THF(RK)*ABE+(.5D0*RK+THG(RK))*ALF1 361 ALF = ALF+2.D0*(.75D0-PA)*THF(RK)*RK*ALF1 362C Viene usato ALF2 per la prima volta 363 ALF = ALF-RK*RK*(PA-.5D0)*ALF2 364C Coefficiente di portanza corretto per Theodorsen (?) 365 CLIFT = ASLRF*ALF 366 OUTA(5) = CLIFT 367 OUTA(6) = CDRAG 368 OUTA(7) = CMOME 369 OUTA(13) = ALF/DEGRAD 370C 371 RETURN 372C 373C 374 END 375C************************* COE2 ********************** 376C= COMPILER (LINK=IBJ$) 377 SUBROUTINE COE2(VCSTR, OUTA, CS, CORDA, RSPEED, 378 & CLIFT, CDRAG, CMOME, ASLOP, BSLOP, CSLOP, 379 & CRF, CFSLOP, DCPDM, DCRDRM, DCMDM, DCRFDM, JPRO) 380 381 implicit none 382 integer*4 i 383C IMPLICIT REAL*8(A-H,O-Z) 384C 385C Input: 386C 387C VCSTR(*) Speed 388C OUTA(*) Array with different parameters, used as I/O 389C CS ? 390C CORDA Chord 391C RSPEED ? 392C 393 real*8 VCSTR(*),CS,CORDA,RSPEED 394 integer*4 JPRO 395C 396C Output: 397C 398C OUTA(*) Array with different parameters, used as I/O 399C CLIFT ? 400C CDRAG ? 401C CMOME ? 402C ASLOP ? 403C BSLOP ? 404C CSLOP ? 405C CRF ? 406C CFSLOP ? 407C DCPDM ? 408C DCRDRM ? 409C DCMDM ? 410C DCRFDM ? 411C DUM ? 412C 413 real*8 OUTA(*),CLIFT,CDRAG,CMOME,ASLOP,BSLOP,CSLOP,CRF,CFSLOP, 414 & DCPDM,DCRDRM,DCMDM,DCRFDM,DUM 415C 416C Local: 417 real*8 ALFA,ALF1,ALF2,VP,VP2,A,B,VC1,GAM,COSGAM,RMACH,ETA,SEGNO, 418 & SEGNOE, 419 & A2,B2,ASN,ASM,SGN,SGM,SGMAX,S2,DAN,DCN,DAM,DCM, 420 & ASLOP0,CSLOP0,ASLRF,C1, 421 & ATMP 422C 423C Data: 424 real*8 PG,DEGRAD,ASN0,ASM0,PN(14),QN(14),PM(14),QM(14) 425 DATA PG /3.14159265358979323846D0/, DEGRAD /.017453293D0/ 426 DATA ASN0 /.22689D0/, ASM0 /.22689D0/ 427 DATA PN /-3.464003D-1,-1.549076D+0, 4.306330D+1,-5.397529D+1, 428 & 5.781402D+0,-3.233003D+1,-2.162257D+1, 1.866347D+1, 429 & 4.198390D+1, 3.295461D+2, 4*0.D0/ 430 DATA QN / 1.533717D+0, 6.977203D+0, 1.749010D+3, 1.694829D 3, 431 & -1.771899D+3,-3.291665D+4, 2.969051D+0,-3.632448D+1, 432 & -2.268578D+3, 6.601995D+3,-9.654208D+3, 8.533930D+4, 433 & -1.492624D+0, 1.163661D+1/ 434 DATA PM / 1.970065D+1,-6.751639D+1, 7.265269D+2, 4.865945D+4, 435 & 2.086279D+4, 6.024672D+3, 1.446334D+2, 8.586896D+2, 436 & -7.550329D+2,-1.021613D+1, 2.247664D+1, 3*0.D0/ 437 DATA QM /-2.322808D+0,-1.322257D+0,-2.633891D+0,-2.180321D-1, 438 & 4.580014D+0, 3.125497D-1,-2.828806D+1,-4.396734D+0, 439 & 2.565870D+2,-1.204976D+1,-1.157802D+2, 8.612138D+0, 440 & 2*0.D0/ 441C 442C CALCOLA I COEFFICIENTI AERODINAMICI CON TEORIA INSTAZIONARIA 443C CON DATI SINTETIZZATI DI BIELAWA 31TH A.H.S. FORUM 1975 444C 445 ALFA = DATAN2(-VCSTR(2),VCSTR(1)) 446 OUTA(2) = ALFA/DEGRAD 447C 448C Questa operazione e' stata spostata all'esterno (COEPRD) per consentire 449C l'uso di integratori impliciti, per i quali la correzione deve essere fatta 450C solo a processo iterativo concluso 451C 452C Coefficienti usati per il ritardo instazionario 453 ALF1 = OUTA(9) 454 ALF2 = OUTA(10) 455C 456 VP2 = VCSTR(1)**2+VCSTR(2)**2 457 VP = DSQRT(VP2) 458 459C print *,'CORDA=',CORDA,', ALF1=',ALF1,', ALF2=',ALF2, 460C & ', VP=',VP,', VP2=',VP2 461 462 A = .5D0*CORDA*ALF1/VP 463 B = .25D0*CORDA*CORDA*ALF2/VP2 464 VC1 = DABS(VCSTR(1)) 465 GAM = DATAN2(-VCSTR(3),VC1) 466 OUTA(3) = GAM/DEGRAD 467 IF(DABS(GAM).GT.PG/3.D0) GAM = PG/3.D0 468 COSGAM = DCOS(GAM) 469 RMACH = DSQRT((VP2+VCSTR(3)**2)*COSGAM)/CS 470 OUTA(4) = RMACH 471 IF(RMACH.GT..99D0) RMACH = .99D0 472 ETA = DSQRT((A/.048D0)**2+(B/.016D0)**2) 473 SEGNO = 1.D0 474 IF(ALFA.LT.0D0) SEGNO = -1.D0 475 SEGNOE = SEGNO 476 IF(ETA.GT.1.D0) SEGNOE = SEGNOE/ETA 477 A = SEGNOE*A 478 B = SEGNOE*B 479 A2 = A*A 480 B2 = B*B 481 ASN = ASN0*(1.D0-RMACH) 482 ASM = ASM0*(1.D0-RMACH) 483 SGN = DABS(ALFA/ASN) 484 SGM = DABS(ALFA/ASM) 485 SGMAX = 1.839D0-70.33D0*B 486 IF(SGMAX.GT.1.86D0) SGMAX = 1.86D0 487 IF(SGN.GT.SGMAX) SGN = SGMAX 488 IF(SGM.GT.SGMAX) SGM = SGMAX 489 490C print *,'A=',A,', B=',B,', SGN=',SGN,', PN(1:6)=',(PN(i),i=1,6) 491 492 DAN = A*(PN(1)+PN(5)*SGN)+B*(PN(2)+PN(6)*SGN) 493 494C print *,'1) DAN=',DAN 495 496 DAN = DAN+DEXP(-1072.52D0*A2)*(A*(PN(3)+PN(7)*SGN)+ 497 & A2*(PN(9)+PN(10)*SGN)) 498 499C print *,'2) DAN=',DAN 500 501 DAN = DAN+DEXP(-40316.42D0*B2)*B*(PN(4)+PN(8)*SGN) 502 503C print *,'3) DAN=',DAN 504 505 DAN = DAN*ASN 506 507C print *,'4) DAN=',DAN 508 509 DCN = A*(QN(1)+QN(3)*A2+SGN*(QN(7)+QN(9)*A2+QN(13)*SGN)+ 510 & B2*(QN(5)+QN(11)*SGN)) 511 DCN = DCN+B*(QN(2)+QN(4)*A2+SGN*(QN(8)+QN(10)*A2+QN(14)*SGN)+ 512 & B2*(QN(6)+QN(12)*SGN)) 513 DAM = A*(PM(1)+PM(3)*A2+PM(5)*B2+PM(10)*SGM+PM(7)*A) 514 DAM = DAM+B*(PM(2)+PM(4)*B2+PM(6)*A2+PM(11)*SGM+PM(8)*B+PM(9)*A) 515 DAM = DAM*ASM 516 S2 = SGM*SGM 517 DCM = A*(QM(2)+QM(8)*A+SGM*(QM(4)+QM(10)*A)+S2*(QM(6)+QM(12)*A)) 518 DCM = DCM+B*(QM(1)+QM(7)*B+SGM*(QM(3)+QM(9)*B)+ 519 & S2*(QM(5)+QM(11)*B)) 520 OUTA(11) = DAN/DEGRAD 521 OUTA(12) = DAM/DEGRAD 522 OUTA(13) = DCN 523 OUTA(14) = DCM 524 525 DAN = SEGNO*DAN 526 DCN = SEGNO*DCN 527 DAM = SEGNO*DAM 528 DCM = SEGNO*DCM 529 530C 531C Questa parte e' sbagliata, perche' per calcolare i coeff. di momento 532C e le loro derivate vengono usati i coeff di portanza e le loro derivate 533C ecc, quindi passando DUM come segnaposto, si ottengono risultati 534C unpredictable. 535C 536C print *,'ATMP = ALFA(=',ALFA,')-DAN(',DAN,')' 537 538 ATMP = ALFA-DAN 539 CALL CPCRCM(ATMP, RMACH, 540 & CLIFT, CDRAG, DUM, ASLOP, BSLOP, DUM, 541 & CRF, CFSLOP, DCPDM, DCRDRM, DUM, DCRFDM, 542 & ASLOP0, CSLOP0, JPRO, 0) 543 544 ATMP = ALFA-DAM 545 CALL CPCRCM(ATMP, RMACH, 546 & DUM, DUM, CMOME, DUM, DUM, CSLOP, 547 & DUM, DUM, DUM, DUM, DCMDM, DUM, 548 & ASLOP0, CSLOP0, JPRO, 1) 549 550 CMOME = CMOME-.25D0*CLIFT 551 CSLOP = CSLOP-.25D0*ASLOP 552 DCMDM = DCMDM-.25D0*DCPDM 553 554 ASLRF = ASLOP0 555 IF(DABS(ALFA-DAN).LT.1.D-6) GOTO 10 556 ASLRF = CLIFT/((ALFA-DAN)*COSGAM) 557 IF(ASLRF.GT.ASLOP0) ASLRF = ASLOP0 558 10 CLIFT = ASLRF*(ALFA-DAN) 559 C1 = .9457D0/DSQRT(1.D0-RMACH*RMACH) 560 561C print *,'COE2 ',CLIFT,ASLOP0*DAN,DCN*C1,CMOME,CSLOP0*DAM,DCM*C1 562 563C La parte steady di CLIFT e' 0 quando alfa cambia segno!!! 564C CLIFT = CLIFT+SEGNO*(ASLOP0*DAN+DCN*C1) 565C CMOME = SEGNO*(CMOME+CSLOP0*DAM+DCM*C1) 566 CLIFT = CLIFT+ASLOP0*DAN+DCN*C1 567 CMOME = CMOME+CSLOP0*DAM+DCM*C1 568 OUTA(5) = CLIFT 569 OUTA(6) = CDRAG 570 OUTA(7) = CMOME 571C 572 RETURN 573C 574C 575 END 576 577 578 579C************************* CPCRCM ******************** 580C= COMPILER (LINK=IBJ$) 581 SUBROUTINE CPCRCM (APHIJ, EMIJ, CLIFT, CDRAG, CMOME, ASLOP, 582 & BSLOP,CSLOP, CRF, CFSLOP, DCPDM, DCDRDM, DCMDM, DCRFDM, 583 & AS0, CS0, JPRO, IRETRN) 584C 585C APHIJ ALFA 586C EMIJ RMACH 587C CLIFT CLIFT 588C CDRAG CDRAG 589C CMOME CMOME 590C ASLOP ASLOP 591C BSLOP BSLOP 592C CSLOP CSLOP 593C CRF CRF 594C CFSLOP CFSLOP 595C DCPDM DCPDM 596C DCDRDM DCDRDM 597C DCMDM DCMDM 598C DCRFDM DCRFDM 599C AS0 ASLOP0 600C CS0 CSLOP0 601C JPRO JPRO 602C IRETRN IRETRN 603C 604 IMPLICIT NONE 605C 606 REAL*8 APHIJ, EMIJ, CLIFT, CDRAG, CMOME, ASLOP, 607 & BSLOP,CSLOP, CRF, CFSLOP, DCPDM, DCDRDM, DCMDM, DCRFDM, 608 & AS0, CS0 609 INTEGER*4 JPRO, IRETRN 610C 611 REAL*8 D1, D2, D3, CD0, B1, B2, B3, 612 & PG, SQT, C1, C2, C3, C4, C5, S, S2, S3, S4, 613 & P1, P2, P3, R0, R1, R2, R3, PEND0, PEND1, PEND2, PEND3, 614 & APHIJ2, APHIJ3 615C 616 INTEGER*4 NEG 617C 618C 619 DIMENSION D1(2,5), D2(2,5), D3(2,5), CD0(2,5), 620 & B1(2,5), B2(2,5), B3(2,5) 621 DATA D1 /.3D0,-.28532D2,.4D0,-.37392D2,.45D0,-.4597D2,.5D0, 622 & -.55819D2,.55D0,-.88102D2/ 623 DATA D2 /.3D0,.4826D1,.4D0,.70982D1,.45D0,.8914D1,.5D0, 624 & .10636D2,.55D0,.1993D2/ 625 DATA D3 /.3D0,.62085D1,.4D0,.64158D1,.45D0,.65334D1,.5D0, 626 & .66613D1,.55D0,.64065D1/ 627 DATA CD0 /.3D0,.60145D-2,.4D0,.35772D-2,.45D0,.2612D-2,.5D0, 628 & .46674D-2,.55D0,.29977D-2/ 629 DATA B1 /.3D0,.86845D-1,.4D0,.1915D0,.45D0,.24923D0,.5D0, 630 & .14568D0,.55D0,.25389D0/ 631 DATA B2 /.3D0,-.91795D0,.4D0,-.1966D1,.45D0,-.26795D1,.5D0, 632 & -.15943D1,.55D0,-.30426D1/ 633 DATA B3 /.3D0,.40357D1,.4D0,.67937D1,.45D0,.91332D1,.5D0, 634 & .62569D1,.55D0,.12049D2/ 635 DATA PG /3.14159265358979323846D0/ 636C**** APHIJ = ANGOLO INCIDENZA 637C**** EMIJ = NUMERO DI MACH 638C**** CLIFT = COEFFICIENTE DI PORTANZA 639C**** CDRAG = COEFFICIENTE DI RESISTENZA TOTALE 640C**** CMOME = COEFFICIENTE DI MOMENTO (CONVENZIONE CABRANTE) RIFERITO 641C**** AL CENTRO AERODINAMICO 642C**** ASLOP = PENDENZA CURVA DI PORTANZA 643C**** BSLOP = PENDENZA CURVA DI RESISTENZA TOTALE 644C**** CSLOP = PENDENZA CURVA DI MOMENTO 645C**** CRF = COEFFICIENTE DI RESISTENZA ATTRITO 646C**** CFSLOP = PENDENZA CURVA DI RESISTENZA ATTRITO 647C**** DCPDM = DERIVATA DI CLIFT RISPETTO NUMERO DI MACH 648C**** DCDRDM = DERIVATA DI CDRAG RISPETTO NUMERO DI MACH 649C**** DCMDM = DERIVATA DI CMOME RISPETTO NUMERO DI MACH 650C**** DCRFDM = DERIVATA DI CRF RISPETTO NUMERO DI MACH 651C**** AS0 = PENDENZA CURVA PORTANZA PER ALFA=0 652C**** CS0 = PENDENZA CURVA MOMENTO PER ALFA=0 653 CFSLOP = 0.D0 654 CRF = .006D0 655 DCRFDM = 0.D0 656 NEG = 1 657 IF(EMIJ.GT..99D0) EMIJ = .99D0 658 SQT = DSQRT(1.D0-EMIJ*EMIJ) 659 C1 = 1.D0-EMIJ 660 C2 = .22689D0*C1 661 C5 = EMIJ/(SQT*SQT) 662 663C print *,'APHIJ: ',APHIJ 664 66597 IF(APHIJ) 181, 182, 182 666181 APHIJ = -APHIJ 667 NEG = -1*NEG 668182 IF(APHIJ-PG) 184, 184, 183 669183 APHIJ = APHIJ-PG*2.D0 670 GOTO 97 671184 IF(APHIJ-C2) 185, 187, 187 672185 ASLOP = 5.7296D0/SQT 673 CLIFT = ASLOP*APHIJ 674 CDRAG = .006D0+.13131D0*APHIJ*APHIJ 675 BSLOP = .26262D0*APHIJ 676 CMOME = 1.4324D0*APHIJ/SQT 677 CSLOP = 1.4324D0/SQT 678 DCPDM = CLIFT*C5 679 DCDRDM = 0.D0 680 DCMDM = CMOME*C5 681 GOTO 250 682187 IF(APHIJ-.34906D0) 189, 191, 191 683189 CLIFT = .29269D0*C1+(1.3D0*EMIJ-.59D0)*APHIJ 684 C2 = (.12217D0+.22689D0*EMIJ)*SQT 685 CMOME = CLIFT/(4*C2) 686 CSLOP = (1.3D0*EMIJ-.59D9)/(4.D0*C2) 687 CLIFT = CLIFT/C2 688 ASLOP = (1.3D0*EMIJ-.59D0)/C2 689 DCPDM = (-.29269D0+1.3D0*APHIJ 690 & -CLIFT*(-(.12217D0+.22689D0*EMIJ)*EMIJ/SQT+ 691 & .22689D0*SQT))/C2 692 DCMDM = .25D0*DCPDM 693 GOTO 210 694 191 IF(APHIJ-2.7402D0) 193, 195, 195 695 193 S = DSIN(APHIJ) 696 S2 = DSIN(2.*APHIJ) 697 S3 = DSIN(3.*APHIJ) 698 S4 = DSIN(4.*APHIJ) 699 CLIFT = (.080373D0*S+1.04308D0*S2 700 & -.011059D0*S3+.023127D0*S4)/SQT 701 CMOME = (-.02827D0*S+.14022D0*S2 702 & -.00622D0*S3+.01012D0*S4)/SQT 703 C1 = DCOS(APHIJ) 704 C2 = DCOS(2.D0*APHIJ) 705 C3 = DCOS(3.D0*APHIJ) 706 C4 = DCOS(4.D0*APHIJ) 707CCCCC CLIFT = CLIFT/C2 708 CSLOP = (-.02827D0*C1+.28044D0*C2 709 & -.01866D0*C3+.04048D0*C4)/SQT 710 ASLOP = (.080373D0*C1+2.08616D0*C2 711 & -.033177D0*C3+.092508D0*C4)/SQT 712 DCPDM = CLIFT*C5 713 DCMDM = CMOME*C5 714 GOTO 210 715195 IF(APHIJ-3.0020D0) 197, 199, 199 716197 CLIFT = -(.4704D0+.10313D0*APHIJ)/SQT 717 ASLOP = -.10313D0/SQT 718 CMOME = -(.4786D0+.02578D0*APHIJ)/SQT 719 CSLOP = .02578D0/SQT 720 DCPDM = CLIFT*C5 721 DCMDM = CMOME*C5 722 GOTO 210 723199 IF(APHIJ-PG) 200, 200, 260 724200 CLIFT = (-17.55D0+5.5864D0*APHIJ)/SQT 725 ASLOP = 5.5864D0/SQT 726 CMOME = (-12.5109D0+3.9824D0*APHIJ)/SQT 727 CSLOP = 3.9824D0/SQT 728 DCPDM = CLIFT*C5 729 DCMDM = CMOME*C5 730210 CDRAG = 1.1233D0-.029894D0*DCOS(APHIJ) 731 & -1.00603D0*DCOS(2.D0*APHIJ) 732 CDRAG = CDRAG+.003115D0*DCOS(3.D0*APHIJ) 733 & -.091487D0*DCOS(4.D0*APHIJ) 734 CDRAG = CDRAG/SQT 735 BSLOP = .029894D0*DSIN(APHIJ)+2.01206D0*DSIN(2.D0*APHIJ) 736 BSLOP = (BSLOP+.009345D0*DSIN(3.D0*APHIJ) 737 & +.365948D0*DSIN(4.D0*APHIJ))/SQT 738 DCDRDM = CDRAG*C5 739 740 250 GOTO(251, 252), JPRO 741C 742C CALCOLO COEFFICIENTI AERODINAMICI PER PROFILO 'RAE9671' 743C 744 252 IF(DABS(APHIJ).GT..21956242D0) GOTO 251 745 IF(NEG.LT.0) APHIJ = -APHIJ 746 CALL LININT(D1, 5, EMIJ, P1, PEND1) 747 CALL LININT(D2, 5, EMIJ, P2, PEND2) 748 CALL LININT(D3, 5, EMIJ, P3, PEND3) 749 APHIJ2 = APHIJ*APHIJ 750 APHIJ3 = APHIJ*APHIJ2 751 CLIFT = P1*APHIJ3+P2*APHIJ2+P3*APHIJ 752 ASLOP = 3.D0*P1*APHIJ2+2.D0*P2*APHIJ+P3 753 AS0 = P3 754 DCPDM = PEND1*APHIJ3+PEND2*APHIJ2+PEND3*APHIJ 755C 756 CALL LININT(CD0, 5, EMIJ, R0, PEND0) 757 CRF = R0 758 CFSLOP = 0.D0 759 DCRFDM = PEND0 760 CALL LININT(B1, 5, EMIJ, R1, PEND1) 761 CALL LININT(B2, 5, EMIJ, R2, PEND2) 762 CALL LININT(B3, 5, EMIJ, R3, PEND3) 763 APHIJ = DABS(APHIJ) 764 APHIJ3 = DABS(APHIJ3) 765 CDRAG = R0+APHIJ*R1+APHIJ2*R2+APHIJ3*R3 766 BSLOP = R1+2.D0*R2*APHIJ+3.D0*R3*APHIJ2 767 DCDRDM = PEND0+APHIJ*PEND1+APHIJ2*PEND2+APHIJ3*PEND3 768C 769 IF(NEG.GT.0) GOTO 249 770 CMOME = -CMOME 771 CSLOP = -CSLOP 772 DCMDM = -DCMDM 773 249 CONTINUE 774 CS0 = 0.D0 775C 776 GOTO 270 777C 778C CMOME = CMOME-.25D0*CLIFT 779C CSLOP = CSLOP-.25D0*ASLOP 780C DCMDM = DCMDM-.25D0*DCPDM 781C RETURN 782C 783C 784C CALCOLO COEFFICIENTI AERODINAMICI PER PROFILO 'NACA0012' 785C 786 251 CONTINUE 787 788 IF(NEG) 255, 255, 260 789255 CLIFT = -CLIFT 790 CMOME = -CMOME 791 APHIJ = -APHIJ 792 DCPDM = -DCPDM 793 DCMDM = -DCMDM 794260 CONTINUE 795C 796 AS0 = 5.7296D0/SQT 797C 798C Qui arriva anche l'altro profilo per una parte di operazioni in comune 799C 800270 CONTINUE 801 CS0 = 0.D0 802 IF(IRETRN.EQ.1) RETURN 803C 804 CMOME = CMOME-.25D0*CLIFT 805 CSLOP = CSLOP-.25D0*ASLOP 806 DCMDM = DCMDM-.25D0*DCPDM 807C 808 RETURN 809 END 810C************************* PSIROT ******************** 811C= COMPILER (LINK=IBJ$) 812 SUBROUTINE PSIROT(PSI, A, NRDA, NC, B, NRDB) 813 IMPLICIT NONE 814 815 INTEGER*4 NRDA, NC, NRDB 816 REAL*8 PSI, A(NRDA,1), B(NRDB,1) 817 818 REAL*8 SINP, COSP, TMP1, TMP2 819 INTEGER*4 K 820 821 822 SINP = DSIN(PSI) 823 COSP = DCOS(PSI) 824 DO 20 K = 1,NC 825 TMP1 = +COSP*A(1,K)-SINP*A(2,K) 826 TMP2 = +SINP*A(1,K)+COSP*A(2,K) 827 B(1,K) = TMP1 828 B(2,K) = TMP2 829 B(3,K) = A(3,K) 830 20 CONTINUE 831 30 CONTINUE 832 RETURN 833 END 834C************************* THF F ********************* 835C= COMPILER (LINK=IBJ$) 836 FUNCTION THF(K) 837 IMPLICIT NONE 838 REAL*8 THF, K, A, B, C, D 839C **** FUNZIONE F DI THEODORSEN *** 840 A = K**4-.102058D0*K**2+9.55732D-6 841 B = -.761036D0*K**3+2.55067D-3*K 842 C = 2.D0*K**4-.113928D0*K**2+9.55732D-6 843 D = -1.064D0*K**3+2.6268D-3*K 844 THF = (A*C+B*D)/(C*C+D*D) 845 RETURN 846 END 847C************************* THG F ********************* 848C= COMPILER (LINK=IBJ$) 849 FUNCTION THG(K) 850 IMPLICIT NONE 851 REAL*8 THG, K, A, B, C, D 852C **** FUNZIONE G DI THEODORSEN *** 853 A = K**4-.102058D0*K**2+9.55732D-6 854 B = -.761036D0*K**3+2.55067D-3*K 855 C = 2.D0*K**4-.113928D0*K**2+9.55732D-6 856 D = -1.064D0*K**3+2.6268D-3*K 857 THG = (B*C-A*D)/(C*C+D*D) 858 RETURN 859 END 860C************************* UNST0 ********************* 861C= COMPILER (LINK=IBJ$) 862 SUBROUTINE UNST0(ALF, PP, COE, PUNTI, NUPIA, NUPIR) 863 IMPLICIT NONE 864 865 INTEGER*4 NUPIA, NUPIR 866 REAL*8 ALF(25,3), PP(5,5), COE(5,5), PUNTI(2,25), 867 & PERM(5), CD(5) 868 869 REAL*8 DSCALL 870 REAL*8 QPSI1, QPSI2,QPSI3 871 INTEGER*4 I, K, L, INDER 872 873 DO 10 I = 1,NUPIA 874 PP(I,1) = 1.D0 875 10 CONTINUE 876 DO 20 K = 2,NUPIA 877 DO 20 I = 1,NUPIA 878 L = (I-1)*NUPIR+1 879 PP(I,K) = PP(I,K-1)*PUNTI(2,L) 880 20 CONTINUE 881 CALL DRUFCT(PP, PERM, NUPIA, 5, INDER) 882 DO 50 K = 1,NUPIR 883 DO 30 I = 1,NUPIA 884 L = (I-1)*NUPIR+K 885 COE(I,1) = ALF(L,1) 886 30 CONTINUE 887 CALL DRUSOL(PP, COE, 5, NUPIA, PERM) 888 DO 40 I = 1,NUPIA 889 L = (I-1)*NUPIR+K 890 QPSI1 = PUNTI(2,L) 891 QPSI2 = QPSI1*QPSI1 892 QPSI3 = QPSI2*QPSI1 893 CD(1) = 0.D0 894 CD(2) = 1.D0 895 CD(3) = 2.D0*QPSI1 896 CD(4) = 3.D0*QPSI2 897 CD(5) = 4.D0*QPSI3 898 ALF(L,2) = DSCALL(CD,COE,NUPIA) 899C 900 CD(1) = 0.D0 901 CD(2) = 0.D0 902 CD(3) = 2.D0 903 CD(4) = 6.D0*QPSI1 904 CD(5) = 12.D0*QPSI2 905 ALF(L,3) = DSCALL(CD,COE,NUPIA) 906 40 CONTINUE 907 50 CONTINUE 908 RETURN 909 END 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933********************************* STLIB ******************************* 934C************************* ADD *********************** 935C= COMPILER (LINK=IBJ$) 936 SUBROUTINE ADD(A,NRDA,B,NRDB,NR,NC,PLUS) 937 IMPLICIT NONE 938 939 INTEGER*4 NRDA, NRDB, NR, NC 940 REAL*8 A(NRDA,1),B(NRDB,1),PLUS 941 942 INTEGER*4 I, K 943 944 DO 10 I=1,NR 945 DO 10 K=1,NC 946 10 A(I,K)=A(I,K)+PLUS*B(I,K) 947 RETURN 948 END 949 950C************************* DABCT3 ******************** 951C= COMPILER (LINK=IBJ$) 952 SUBROUTINE DABCT3(A,B,NRDB,C) 953C 954C CALCOLA IL PRODOTTO DELLE MATRICI (3*3) "A" * "B" * "C-TRASPOSTA" 955C 956 IMPLICIT NONE 957 958 INTEGER*4 NRDB 959 REAL*8 A,B,C,D 960 961 DIMENSION A(3,3),B(NRDB,1),C(3,3),D(3,3) 962 963 INTEGER*4 I, K, L 964 965 DO 10 I=1,3 966 DO 10 K=1,3 967 D(I,K)=0.D0 968 DO 10 L=1,3 96910 D(I,K)=D(I,K)+A(I,L)*B(L,K) 970 DO 20 I=1,3 971 DO 20 K=1,3 972 B(I,K)=0.D0 973 DO 20 L=1,3 97420 B(I,K)=B(I,K)+D(I,L)*C(K,L) 975 RETURN 976 END 977 978C ****************** DDEMA *********************** 979 SUBROUTINE DDEMA(RO,RD,DEMA,ND,FACT) 980C 981C GIVEN THE ROTATION VECTOR RO AND THE DERIVATIVE OF THE ROTATION 982C VECTOR RD, CALCULATE THE CAPITAL GAMMA DOT. 983C 984 IMPLICIT NONE 985 986 INTEGER*4 ND 987 REAL*8 RO,RD,DEMA,FACT 988 989 REAL*8 A,B,C,D,E,G,RDRO,ARG,ARG2,T 990 INTEGER*4 I 991 992 DIMENSION RO(3),RD(3),DEMA(ND,3),T(3) 993 994 REAL*8 TRESH2 995 DATA TRESH2 /1.D-10/ 996 997 ARG2=0.D0 998 RDRO=0.D0 999 DO 10 I=1,3 1000 T(I)=RO(I)*FACT 1001 ARG2=ARG2+T(I)**2 1002 RDRO=RDRO+T(I)*RD(I) 1003 10 CONTINUE 1004 IF (ARG2.LT.TRESH2) THEN 1005C 1006 A=(1.D0-ARG2/30.D0*(2.D0+ARG2/56.D0*(3.D0-ARG2/22.5D0)))/12.D0 1007 B=(1.D0-ARG2/42.D0*(2.D0-ARG2/72.D0*(3.D0-ARG2/27.5D0)))/60.D0 1008 D=(1.D0-ARG2/12.D0*(1.D0-ARG2/30.D0*(1.D0-ARG2/56.0D0)))/ 2.D0 1009 E=(1.D0-ARG2/20.D0*(1.D0-ARG2/42.D0*(1.D0-ARG2/72.0D0)))/ 6.D0 1010 A=-A*RDRO 1011 B=-B*RDRO 1012C 1013 ELSE 1014C 1015 ARG=DSQRT(ARG2) 1016 D=(1.D0-DCOS(ARG))/ARG2 1017 G=DSIN(ARG)/ARG 1018 E=(1.D0-G)/ARG2 1019 A=(G-2.D0*D)*RDRO/ARG2 1020 B=(D-3.D0*E)*RDRO/ARG2 1021C 1022 END IF 1023C 1024 C=-B*ARG2 1025 DEMA(1,1)=T(1)*T(1)*B+C 1026 DEMA(1,2)=T(1)*T(2)*B-T(3)*A 1027 DEMA(1,3)=T(1)*T(3)*B+T(2)*A 1028 DEMA(2,1)=T(2)*T(1)*B+T(3)*A 1029 DEMA(2,2)=T(2)*T(2)*B+C 1030 DEMA(2,3)=T(2)*T(3)*B-T(1)*A 1031 DEMA(3,1)=T(3)*T(1)*B-T(2)*A 1032 DEMA(3,2)=T(3)*T(2)*B+T(1)*A 1033 DEMA(3,3)=T(3)*T(3)*B+C 1034 DEMA(1,1)=DEMA(1,1)+ 2.D0*(RD(1)*T(1)-RDRO)*E 1035 DEMA(1,2)=DEMA(1,2)+(RD(1)*T(2)+T(1)*RD(2))*G-RD(3)*D 1036 DEMA(1,3)=DEMA(1,3)+(RD(1)*T(3)+T(1)*RD(3))*E+RD(2)*D 1037 DEMA(2,1)=DEMA(2,1)+(RD(2)*T(1)+T(2)*RD(1))*E+RD(3)*D 1038 DEMA(2,2)=DEMA(2,2)+ 2.D0*(RD(2)*T(2)-RDRO)*E 1039 DEMA(2,3)=DEMA(2,3)+(RD(2)*T(3)+T(2)*RD(3))*E-RD(1)*D 1040 DEMA(3,1)=DEMA(3,1)+(RD(3)*T(1)+T(3)*RD(1))*E-RD(2)*D 1041 DEMA(3,2)=DEMA(3,2)+(RD(3)*T(2)+T(3)*RD(2))*E+RD(1)*D 1042 DEMA(3,3)=DEMA(3,3)+ 2.D0*(RD(3)*T(3)-RDRO)*E 1043 RETURN 1044 END 1045 1046C ****************** DLGMA *********************** 1047 SUBROUTINE DLGMA(RO,SPMA,ND,FACT,ID) 1048C 1049C GIVEN THE ROTATION VECTOR RO, CALCULATE THE CAPITAL GAMMA MATRIX 1050C SPMA. 1051C 1052 IMPLICIT NONE 1053 1054 INTEGER*4 ND, ID 1055 REAL *8 RO,SPMA,FACT 1056 1057 REAL*8 A,B,C,ARG,ARG2,T,TRESH2 1058 INTEGER*4 I 1059 1060 DIMENSION RO(3),SPMA(ND,3),T(3) 1061 1062 DATA TRESH2 /1.D-10/ 1063 1064 ARG2=0.D0 1065 DO 10 I=1,3 1066 T(I)=RO(I)*FACT 1067 ARG2=ARG2+T(I)**2 1068 10 CONTINUE 1069 IF (ARG2.LT.TRESH2) THEN 1070 A=(1.D0-ARG2/12.D0*(1-ARG2/30.D0*(1-ARG2/56.D0)))/2.D0 1071 B=(1.D0-ARG2/20.D0*(1-ARG2/42.D0*(1-ARG2/72.D0)))/6.D0 1072 ELSE 1073 ARG=DSQRT(ARG2) 1074 A=(1.D0-DCOS(ARG))/ARG2 1075 B=(1.D0-DSIN(ARG)/ARG)/ARG2 1076 END IF 1077 IF (ID.LT.0) THEN 1078 IF (ARG2.LT.TRESH2) THEN 1079 B=(1.D0+ARG2/60.D0*(1.D0+ARG2/420.D0))/12.D0 1080 ELSE 1081 B=(1.D0-.5D0*(1.D0-B*ARG2)/A)/ARG2 1082 END IF 1083 A=-.5D0 1084 END IF 1085 C=1.D0-B*ARG2 1086 SPMA(1,1)=T(1)*T(1)*B+C 1087 SPMA(1,2)=T(1)*T(2)*B-T(3)*A 1088 SPMA(1,3)=T(1)*T(3)*B+T(2)*A 1089 SPMA(2,1)=T(2)*T(1)*B+T(3)*A 1090 SPMA(2,2)=T(2)*T(2)*B+C 1091 SPMA(2,3)=T(2)*T(3)*B-T(1)*A 1092 SPMA(3,1)=T(3)*T(1)*B-T(2)*A 1093 SPMA(3,2)=T(3)*T(2)*B+T(1)*A 1094 SPMA(3,3)=T(3)*T(3)*B+C 1095 RETURN 1096 END 1097 1098C************************* DPRMMA ******************** 1099C= COMPILER (LINK=IBJ$) 1100 SUBROUTINE DPRMMA(A,NRDA,NRA,NCA,B,NRDB,C,NRDC,NCC,IA,IB) 1101C 1102C ESEGUE IL PRODOTTO DELLA MATRICE "A" PER "B" E SOMMA 1103C IL RISULTATO IN "C" 1104C A = PRIMA MATRICE 1105C NRDA = N.RO DI RIGHE PER IL DIMENSIONAMENTO DI "A" 1106C NRA = N.RO DI RIGHE UTILIZZATE PER IL PRODOTTO DI "A" 1107C NCA = N.RO DI COLONNE DI "A" 1108C B = SECONDA MATRICE 1109C NRDB = N.RO DI RIGHE PER IL DIMENSIONAMENTO DI "B" 1110C C = MATRICE PRODOTTO 1111C NRDC = N.RO DI RIGHE PER IL DIMENSIONAMENTO DI "C" 1112C NCB = N.RO DI COLONNE DI "C" 1113C IA = INDICE DI TRASPOSIZIONE DI "A" 0 = NO 1114C 1 = SI 1115C IB = INDICE DI TRASPOSIZIONE DI "B" 1116C 1117 IMPLICIT NONE 1118 1119 INTEGER*4 NRDA,NRA,NCA,NRDB,NRDC,NCC,IA,IB 1120 REAL*8 A,B,C 1121 DIMENSION A(NRDA,1),B(NRDB,1),C(NRDC,1) 1122 1123 INTEGER*4 IGO,I,K,L 1124 1125 IGO=2*IA+IB+1 1126 GOTO(10,20,30,40),IGO 112710 DO 11 I=1,NRA 1128 DO 11 K=1,NCC 1129 DO 11 L=1,NCA 1130 C(I,K)=C(I,K)+A(I,L)*B(L,K) 113111 CONTINUE 1132 RETURN 113320 DO 12 I=1,NRA 1134 DO 12 K=1,NCC 1135 DO 12 L=1,NCA 1136 C(I,K)=C(I,K)+A(I,L)*B(K,L) 113712 CONTINUE 1138 RETURN 113930 DO 13 I=1,NCA 1140 DO 13 K=1,NCC 1141 DO 13 L=1,NRA 1142 C(I,K)=C(I,K)+A(L,I)*B(L,K) 114313 CONTINUE 1144 RETURN 114540 DO 14 I=1,NCA 1146 DO 14 K=1,NCC 1147 DO 14 L=1,NRA 1148 C(I,K)=C(I,K)+A(L,I)*B(K,L) 114914 CONTINUE 1150 RETURN 1151 END 1152 1153C************************* DPROMM ******************** 1154C= COMPILER (LINK=IBJ$) 1155 SUBROUTINE DPROMM(A,NRDA,NRA,NCA,B,NRDB,C,NRDC,NCC,IA,IB) 1156C 1157C ESEGUE IL PRODOTTO DELLA MATRICE "A" PER "B" 1158C A = PRIMA MATRICE 1159C NRDA = N.RO DI RIGHE PER IL DIMENSIONAMENTO DI "A" 1160C NRA = N.RO DI RIGHE UTILIZZATE PER IL PRODOTTO DI "A" 1161C NCA = N.RO DI COLONNE DI "A" 1162C B = SECONDA MATRICE 1163C NRDB = N.RO DI RIGHE PER IL DIMENSIONAMENTO DI "B" 1164C C = MATRICE PRODOTTO 1165C NRDC = N.RO DI RIGHE PER IL DIMENSIONAMENTO DI "C" 1166C NCB = N.RO DI COLONNE DI "C" 1167C IA = INDICE DI TRASPOSIZIONE DI "A" 0 = NO 1168C 1 = SI 1169C IB = INDICE DI TRASPOSIZIONE DI "B" 1170C 1171 IMPLICIT NONE 1172 1173 INTEGER*4 NRDA,NRA,NCA,NRDB,NRDC,NCC,IA,IB 1174 REAL*8 A,B,C 1175 1176 DIMENSION A(NRDA,1),B(NRDB,1),C(NRDC,1) 1177 1178 INTEGER*4 IGO,I,K,L 1179 1180 IGO=2*IA+IB+1 1181 GOTO(10,20,30,40),IGO 118210 DO 11 I=1,NRA 1183 DO 11 K=1,NCC 1184 C(I,K)=0.D0 1185 DO 11 L=1,NCA 1186 C(I,K)=C(I,K)+A(I,L)*B(L,K) 118711 CONTINUE 1188 RETURN 118920 DO 12 I=1,NRA 1190 DO 12 K=1,NCC 1191 C(I,K)=0.D0 1192 DO 12 L=1,NCA 1193 C(I,K)=C(I,K)+A(I,L)*B(K,L) 119412 CONTINUE 1195 RETURN 119630 DO 13 I=1,NCA 1197 DO 13 K=1,NCC 1198 C(I,K)=0.D0 1199 DO 13 L=1,NRA 1200 C(I,K)=C(I,K)+A(L,I)*B(L,K) 120113 CONTINUE 1202 RETURN 120340 DO 14 I=1,NCA 1204 DO 14 K=1,NCC 1205 C(I,K)=0.D0 1206 DO 14 L=1,NRA 1207 C(I,K)=C(I,K)+A(L,I)*B(K,L) 120814 CONTINUE 1209 RETURN 1210 END 1211 1212 1213C************************* DPROMV ******************** 1214C= COMPILER (LINK=IBJ$) 1215 SUBROUTINE DPROMV(A,NRDA,NRA,NCA,B,C,IA) 1216C 1217C ESEGUE IL PRODOTTO DELLA MATRICE "A" PER IL VETTORE "B" 1218C 1219 IMPLICIT NONE 1220 1221 INTEGER*4 NRDA,NRA,NCA,IA 1222 REAL*8 A,B,C 1223 1224 INTEGER*4 IGO,I,K 1225 REAL*8 D 1226 1227 DIMENSION A(NRDA,1),B(NCA),C(NRA) 1228 1229 IGO=IA+1 1230 GOTO(100,200),IGO 1231100 DO 20 I=1,NRA 1232 D=0.D0 1233 DO 10 K=1,NCA 123410 D=D+A(I,K)*B(K) 123520 C(I)=D 1236 RETURN 1237200 DO 40 I=1,NRA 1238 D=0.D0 1239 DO 30 K=1,NCA 124030 D=D+A(K,I)*B(K) 124140 C(I)=D 1242 RETURN 1243 END 1244 1245C************************* DPRVAD ******************** 1246C= COMPILER (LINK=IBJ$) 1247 SUBROUTINE DPRVAD(X,Y,Z) 1248C 1249C ESEGUE IL PRODOTTO DEL VETTORE "X" 1250C PER IL VETTORE "Y" E SOMMA IL RISULTATO IN "Z" 1251C 1252 IMPLICIT NONE 1253 1254 REAL*8 X,Y,Z,W 1255 1256 DIMENSION X(3),Y(3),Z(3),W(3) 1257 1258 W(1)= X(2)*Y(3)-X(3)*Y(2) 1259 W(2)= X(3)*Y(1)-X(1)*Y(3) 1260 W(3)= X(1)*Y(2)-X(2)*Y(1) 1261 Z(1)=Z(1)+W(1) 1262 Z(2)=Z(2)+W(2) 1263 Z(3)=Z(3)+W(3) 1264 RETURN 1265 END 1266 1267C************************* DPRVVA ******************** 1268C= COMPILER (LINK=IBJ$) 1269 SUBROUTINE DPRVVA(A,B,C,NRDC,NRA,NRB) 1270 1271 IMPLICIT NONE 1272 1273 INTEGER*4 NRDC,NRA,NRB 1274 REAL*8 A(1),B(1),C(NRDC,1) 1275 1276 INTEGER*4 I,K 1277C 1278C ESEGUE IL PRODOTTO DEL VETTORE A PER IL VETTORE B TRASPOSTO 1279C LA MATRICE RISULTANTE VIENE SOMMATA IN C 1280C 1281 DO 20 I=1,NRA 1282 DO 20 K=1,NRB 1283 C(I,K)=C(I,K)+A(I)*B(K) 1284 20 CONTINUE 1285C 1286 RETURN 1287C 1288C 1289 END 1290 1291C************************* DRMOVE ******************** 1292C= COMPILER (LINK=IBJ$) 1293 SUBROUTINE DRMOVE(A,B,N,H) 1294C 1295C SOMMA AI TERMINI DEL VETTORE "A" QUELLI DI "B" MOLTIPLICATI PER H 1296C 1297 IMPLICIT NONE 1298 1299 INTEGER*4 N 1300 REAL*8 A,B,H 1301 DIMENSION A(1),B(1) 1302 1303 INTEGER*4 I 1304 1305 DO 10 I=1,N 1306 A(I)=A(I)+B(I)*H 130710 CONTINUE 1308 RETURN 1309 END 1310 1311C************************* DRUFCT ******************** 1312C= COMPILER (LINK=IBJ$) 1313 SUBROUTINE DRUFCT(A,PERM,N,NR,INDER) 1314 IMPLICIT NONE 1315 1316 INTEGER*4 N,NR,INDER 1317 REAL*8 A,PERM 1318 1319 DIMENSION A(NR,1),PERM(1) 1320 1321 REAL*8 X,DP 1322 INTEGER*4 I,K,IM1,IP1,IPVT,J 1323 1324 DO10I=1,N 1325 X=0.D0 1326 DO20K=1,N 1327 IF(DABS(A(I,K)).LT.X)GOTO20 1328 X=DABS(A(I,K)) 132920 CONTINUE 1330 IF(X.EQ.0.D0)GOTO110 1331 PERM(I)=1.D0/X 133210 CONTINUE 1333 DO100 I=1,N 1334 IM1=I-1 1335 IP1=I+1 1336 IPVT=I 1337 X=0.D0 1338 DO50K=I,N 1339 DP=A(K,I) 1340 IF(I.EQ.1)GOTO40 1341 DO30J=1,IM1 1342 DP=DP-A(K,J)*A(J,I) 134330 CONTINUE 1344 A(K,I)=DP 134540 IF(X.GT.(DABS(DP)*PERM(K)))GOTO50 1346 IPVT=K 1347 X=DABS(DP)*PERM(K) 134850 CONTINUE 1349 IF(X.LE.0.D0)GOTO110 1350 IF(IPVT.EQ.I)GOTO70 1351 DO60J=1,N 1352 X=A(IPVT,J) 1353 A(IPVT,J)=A(I,J) 1354 A(I,J)=X 135560 CONTINUE 1356 PERM(IPVT)=PERM(I) 135770 PERM(I)=IPVT 1358 IF(I.EQ.N)GOTO100 1359 X=A(I,I) 1360 DO90K=IP1,N 1361 A(K,I)=A(K,I)/X 1362 IF(I.EQ.1)GOTO90 1363 DP=A(I,K) 1364 DO80J=1,IM1 1365 DP=DP-A(I,J)*A(J,K) 136680 CONTINUE 1367 A(I,K)=DP 136890 CONTINUE 1369100 CONTINUE 1370 INDER=0 1371 RETURN 1372110 INDER=I 1373 RETURN 1374 END 1375 1376C************************* DRUSOL ******************** 1377C= COMPILER (LINK=IBJ$) 1378 SUBROUTINE DRUSOL(A,B,NR,N,PERM) 1379 1380 IMPLICIT NONE 1381 1382 INTEGER*4 NR,N 1383 REAL*8 A,B,PERM 1384 1385 DIMENSION A(NR,1),B(1),PERM(1) 1386 1387 REAL*8 X,DP 1388 INTEGER*4 I,K,IM1,INF 1389 1390 IF(N.GT.1)GOTO 2 1391 B(1)=B(1)/A(1,1) 1392 RETURN 13932 DO10I=1,N 1394 K=PERM(I) 1395 IF(K.EQ.I)GOTO10 1396 X=B(K) 1397 B(K)=B(I) 1398 B(I)=X 139910 CONTINUE 1400 DO20I=2,N 1401 IM1=I-1 1402 DP=B(I) 1403 DO40K=1,IM1 1404 DP=DP-A(I,K)*B(K) 140540 CONTINUE 1406 B(I)=DP 140720 CONTINUE 1408 B(N)=B(N)/A(N,N) 1409 DO60I=2,N 1410 IM1=N-I+1 1411 INF=IM1+1 1412 DP=B(IM1) 1413 DO70K=INF,N 1414 DP=DP-A(IM1,K)*B(K) 141570 CONTINUE 1416 B(IM1)=DP/A(IM1,IM1) 141760 CONTINUE 1418 RETURN 1419 END 1420 1421 1422C************************* DSCALL ********************* 1423C= COMPILER (LINK=IBJ$) 1424 REAL*8 FUNCTION DSCALL(X,Y,N) 1425C 1426C ESEGUE IL PRODOTTO SCALARE DI X PER Y 1427C 1428 IMPLICIT NONE 1429 INTEGER*4 N 1430 REAL*8 X, Y 1431 1432 INTEGER*4 I 1433 DIMENSION X(1),Y(1) 1434 DSCALL=0.D0 1435 DO 1 I=1,N 14361 DSCALL=DSCALL+X(I)*Y(I) 1437 RETURN 1438 END 1439 1440C************************* DVET ********************** 1441C= COMPILER (LINK=IBJ$) 1442 SUBROUTINE DVET(W,WV,IVER) 1443C 1444C CALCOLA L'OPERATORE VETTORE DI UN VETTORE 1445C 1446 IMPLICIT NONE 1447 1448 INTEGER*4 IVER 1449 REAL*8 W,WV 1450 1451 REAL*8 H 1452 1453 DIMENSION W(3),WV(9) 1454 1455 H=DFLOAT(IVER) 1456 WV(1)= 0.D0 1457 WV(2)= H*W(3) 1458 WV(3)=-H*W(2) 1459 WV(4)= -WV(2) 1460 WV(5)= 0.D0 1461 WV(6)= H*W(1) 1462 WV(7)= -WV(3) 1463 WV(8)= -WV(6) 1464 WV(9)= 0.D0 1465 RETURN 1466 END 1467 1468 1469C************************* LININT ******************** 1470C= COMPILER (LINK=IBJ$) 1471 SUBROUTINE LININT(XY,N,X,Y,PEND) 1472 1473 IMPLICIT NONE 1474 1475 INTEGER*4 N 1476 REAL*8 XY(2,1),X,Y,PEND 1477 1478 INTEGER*4 K1, K2, I 1479 IF(X.GT.XY(1,1)) GO TO 10 1480 K1=1 1481 GO TO 30 1482 10 DO 20 I=2,N 1483 K1=I-1 1484 IF(XY(1,I).GT.X) GO TO 30 1485 20 CONTINUE 1486 K1=N-1 1487 30 K2=K1+1 1488 PEND=(XY(2,K2)-XY(2,K1))/(XY(1,K2)-XY(1,K1)) 1489 Y=XY(2,K1)+PEND*(X-XY(1,K1)) 1490 RETURN 1491 END 1492 1493C************************* MOVE ********************** 1494C= COMPILER (LINK=IBJ$) 1495 SUBROUTINE MOVE(A,B,N) 1496C 1497C TRASFERISCE N TERMINI DEL VETTORE 'B' NEL VETTORE 'A' 1498C 1499 IMPLICIT NONE 1500 1501 INTEGER*4 N 1502 REAL*8 A,B 1503 DIMENSION A(1),B(1) 1504 1505 INTEGER*4 I 1506 1507 DO 10 I=1,N 1508 A(I)=B(I) 150910 CONTINUE 1510 RETURN 1511 END 1512 1513 1514C************************* PRM3 ******************** 1515 SUBROUTINE PRM3 (A,NRDA,B,NRDB,C,NRDC,IA,IB) 1516C 1517C ESEGUE IL PRODOTTO DELLA MATRICE "A" PER "B" 1518C IL PRODOTTO E' DI 3 RIGHE PER 3 COLONNE 1519C A = PRIMA MATRICE 1520C NRDA = N.RO DI RIGHE PER IL DIMENSIONAMENTO DI "A" 1521C NCA = N.RO DI COLONNE DI "A" 1522C B = SECONDA MATRICE 1523C NRDB = N.RO DI RIGHE PER IL DIMENSIONAMENTO DI "B" 1524C C = MATRICE PRODOTTO 1525C NRDC = N.RO DI RIGHE PER IL DIMENSIONAMENTO DI "C" 1526C IA = INDICE DI TRASPOSIZIONE DI "A" 0 = NO 1527C 1 = SI 1528C IB = INDICE DI TRASPOSIZIONE DI "B" 1529C 1530 IMPLICIT NONE 1531 1532 INTEGER*4 NRDA,NRDB,NRDC,IA,IB 1533 REAL*8 A,B,C 1534 1535 INTEGER*4 NRA, NCA, NCC, IGO, I, K, L 1536 REAL*8 D 1537 DIMENSION A(NRDA,1),B(NRDB,1),C(NRDC,1),D(3,3) 1538 1539 NRA=3 1540 NCA=3 1541 NCC=3 1542 IGO=2*IA+IB+1 1543 GOTO(10,20,30,40),IGO 154410 DO 11 I=1,NRA 1545 DO 11 K=1,NCC 1546 D(I,K)=0.D0 1547 DO 11 L=1,NCA 1548 D(I,K)=D(I,K)+A(I,L)*B(L,K) 154911 CONTINUE 1550 GO TO 50 155120 DO 12 I=1,NRA 1552 DO 12 K=1,NCC 1553 D(I,K)=0.D0 1554 DO 12 L=1,NCA 1555 D(I,K)=D(I,K)+A(I,L)*B(K,L) 155612 CONTINUE 1557 GO TO 50 155830 DO 13 I=1,NCA 1559 DO 13 K=1,NCC 1560 D(I,K)=0.D0 1561 DO 13 L=1,NRA 1562 D(I,K)=D(I,K)+A(L,I)*B(L,K) 156313 CONTINUE 1564 GO TO 50 156540 DO 14 I=1,NCA 1566 DO 14 K=1,NCC 1567 D(I,K)=0.D0 1568 DO 14 L=1,NRA 1569 D(I,K)=D(I,K)+A(L,I)*B(K,L) 157014 CONTINUE 1571 50 DO 60 I=1,3 1572 DO 60 K=1,3 1573 C(I,K)=D(I,K) 1574 60 CONTINUE 1575 RETURN 1576 END 1577 1578 1579 1580C************************* DZERO ******************** 1581 SUBROUTINE DZERO(A,N) 1582 1583 IMPLICIT NONE 1584 1585 INTEGER*4 N,I 1586 REAL*8 A(N) 1587 1588 DO I=1,N 1589 A(I)=0.D0 1590 END DO 1591 RETURN 1592 END 1593 1594C************************* POLCOE ******************** 1595 SUBROUTINE POLCOE(X,Y,N,COF) 1596 1597 IMPLICIT NONE 1598 1599 REAL*8 X(1),Y(1),COF(1) 1600 INTEGER*4 N 1601C 1602 REAL*8 PHI,B,FF,S(10) 1603 INTEGER*4 I,J,K 1604C 1605C Fixme: check for N <= 10? 1606 DO I=1,N 1607 S(I)=0.D0 1608 COF(I)=0.D0 1609 ENDDO 1610 S(N)=-X(1) 1611 DO I=2,N 1612 DO J=N+1-I,N-1 1613 S(J)=S(J)-X(I)*S(J+1) 1614 ENDDO 1615 S(N)=S(N)-X(I) 1616 ENDDO 1617 DO J=1,N 1618 PHI=N 1619 DO K=N-1,1,-1 1620 PHI=K*S(K+1)+X(J)*PHI 1621 ENDDO 1622 FF=Y(J)/PHI 1623 B=1.D0 1624 DO K=N,1,-1 1625 COF(K)=COF(K)+B*FF 1626 B=S(K)+X(J)*B 1627 ENDDO 1628 ENDDO 1629 RETURN 1630 END 1631 1632