+ 00000010
+ PROGRAM FRANKIE ! ABBIAMO CREATO UN MOSTRO? 00000020
+C ************************************************************* 00000030
+C 00000040
+C QUESTO PROGRAMMA CALCOLA LE POSIZIONI DEI PIANETI E DELLA LUNA 00000050
+C 00000060
+C COMITATO DI REDAZIONE: 00000070
+C MARCELLO GALLI (GRUPPO 5) 00000080
+C ENRICA BAIADA (GRUPPO 5) 00000090
+C SABINA MENGOLI (GRUPPO 8) 00000100
+C THIS ISSUE: SEPT-OCT 84 LAST REVISION:DEC 84 00000110
+C 00000120
+C SEMBRA UN ESERCIZIO DI MECCANICA CELESTE, MA IN REALTA' IL 00000130
+C SUO UNICO BIECO SCOPO E' QUELLO DI PERMETTERCI DI FARE OROSCOPI,00000140
+C 00000150
+C Non e' colpa nostra se il programma e' venuto lungo e 00000160
+C COMPLICATO! LA TERRA SI AGITA TROPPO PER ESSERE UN COMODO 00000170
+C punto di osservazione. 00000180
+C ************************************************************* 00000190
+C 00000200
+ IMPLICIT REAL*8 (A-H,O-Z) 00000210
+ CHARACTER*30,NOME 00000220
+ INTEGER ANNO 00000230
+ INTEGER CASE(10),SEGNI(10) 00000240
+ REAL*8 MS,ML,MT,LAMBDA 00000250
+ COMMON/MASS/MS,MT,ML,AMS(9),G,PIGREC,UA,EPS,PASSOP, 00000260
+ 1 STMP,KRETT 00000270
+ DIMENSION X(3,10),F(9),AE(9),XL(3) 00000280
+ DIMENSION AME(9) ! ANOMALIA MEDIA 00000290
+ DIMENSION BETA(10),LAMBDA(10),ISEGNC(12) 00000300
+ DIMENSION CUSP(12),X1(3,10),ALFA(10),DELTA(10) 00000310
+ DIMENSION E(6,9) 00000320
+ DIMENSION AM(9) 00000330
+C 00000340
+C E sono gli elementi orbitali presi dal file CIBO.DAT preparato 00000350
+C da IGOR.FOR 00000360
+C Gli elementi orbitali usati sono: (la terra e' al primo posto) 00000370
+C inclinazione orbita sull'eclittica 00000380
+C longitudine nodo ascendente (omega piccolo) 00000390
+C longitudine perielio (omega grande) 00000400
+C eccetricita' 00000410
+C Anomalia media all'epoca misurata dal perielio 00000420
+C SEMIASSE MAGGIORE ( A ),IN KM!. ANGOLI CONTATI ORARI DA GAMMA 00000430
+C 00000440
+C AM sono le masse dei pianeti, 00000450
+C MS=massa luna,MS=massa sole,MT=massa terra(le masse sono date in00000460
+C X= posizioni dei pianeti (eliocentriche) 00000470
+C F= anomalia vera dei pianeti 00000480
+C AE=anomalia eccentrica dei pianeti 00000490
+C XL=posizione della luna 00000500
+C CUSP=longituduine eclittica cuspidi delle case 00000510
+C CASE=numero casa dei pianeti 00000520
+C SEGNI=costellazione del pianeta 00000530
+C ISEGNC= segni delle cuspidi delle case 00000540
+C BETA= latitudine eclittica BETL per la luna 00000550
+C LAMBDA= longitudine eclittica ALAMBL per la luna 00000560
+C ALFA=ascensione retta 00000570
+C DELTA=declinazione 00000580
+C X1=coord rettangolari geocentriche pianeti 00000590
+C -------------------------------------------------------- 00000600
+ DATA AM/6.0469D24,3.3019D23,4.8693D24, 00000610
+ 1 6.4248D23,1.899174D27,5.686399D26, 00000620
+ 2 8.6633D25,1.0299D26,1.5913D22/ 00000630
+ PIGREC=3.14159265358979323846 00000640
+C costante gravitazionale IN 1000*NEWTON*Km**2/Kg**2 00000650
+ G=6.6732D-20 00000660
+ UA=149.6D6 00000670
+ MS=1.9891D30 00000680
+ ML=7.36D22 00000690
+ MT=5.9734D24 00000700
+ DO 99 I=1,9 00000710
+ 99 AMS(I)=AM(I) 00000720
+C 00000730
+C 00000740
+C ------------------------------------------------------------ 00000750
+ PASSOP=0.D0 00000760
+ STMP=0.D0 00000770
+ KRETT=0 00000780
+C -------------------- input da terminale : 00000790
+C 00000800
+ 1 CONTINUE 00000810
+ TYPE*,' Come ti chiami? (metti il nome fra apici)' 00000820
+ ACCEPT*,NOME 00000830
+ IF(NOME.NE.'TABELLE') GO TO100 00000840
+ TYPE*,' Fornire i parametri di controllo:' 00000850
+ TYPE*,' STMP,PASSOP,KRETT' 00000860
+ ACCEPT*,STMP,PASSOP,KRETT 00000870
+ TYPE*,' LETTI PARAMETRI DI CONTROLLO:',STMP,PASSOP 00000880
+ 1 ,KRETT 00000890
+ GOTO 1 00000900
+ 100 CONTINUE 00000910
+ TYPE*,' Quando sei nato? (metti giorno mese anno )' 00000920
+ ACCEPT*,GIORNO,MESE,ANNO 00000930
+ IF(TESTG(ANNO,MESE,GIORNO).LT.1.D0) GO TO 2 00000940
+ TYPE*,' MI HAI DATO UN GIORNO CHE NON ESISTE!' 00000950
+ GO TO 1 00000960
+ 2 TYPE*,' A che ora? (metti l''ora nella forma: ore.minuti )' 00000970
+ TYPE*,' TEMPO UNIVERSALE!'
+ TYPE*,' (Che e'' quello di Greenwich! BADA A FUSI E ORA LEGALE)'
+ ACCEPT*,ORA 00000990
+ IF(ORA.LE.24.D0) GO TO 101 00001000
+ TYPE*,' COME ?' 00001010
+ GO TO 2 00001020
+ 101 CONTINUE 00001030
+ TYPE*,' A che latitudine e longitudine? '
+ TYPE*,' LATITUDINE IN GRADI.PRIMI '
+ TYPE*,' LONGITUDINE CONTATA VERSO EST! IN ORE.MINUTI' 00001050
+ ACCEPT*,ALAT,ALONG 00001060
+ TYPE*,'SEMBRA CHE TU SIA:',NOME,' NATO IL:',ANNO,MESE,GIORNO 00001070
+ TYPE*,' ALLE ORE:',ORA,' A:',ALAT,ALONG 00001080
+C Intestazioni dei files di output 00001090
+ IF(STMP.LE.0.D0) GO TO 102 00001100
+ WRITE(2,1000)NOME,ANNO,MESE,GIORNO 00001110
+ WRITE(3,1000)NOME,ANNO,MESE,GIORNO 00001120
+ WRITE(4,1000)NOME,ANNO,MESE,GIORNO 00001130
+ WRITE(11,1000)NOME,ANNO,MESE,GIORNO 00001140
+ WRITE(2,2000)ORA,ALAT,ALONG,STMP,PASSOP,KRETT 00001150
+ WRITE(3,2000)ORA,ALAT,ALONG,STMP,PASSOP,KRETT 00001160
+ WRITE(4,2000)ORA,ALAT,ALONG,STMP,PASSOP,KRETT 00001170
+ WRITE(11,2000)ORA,ALAT,ALONG,STMP,PASSOP,KRETT 00001180
+ 1000 FORMAT(' NOME:',A30,' DATA:',I5,I5,F10.5) 00001190
+ 2000 FORMAT(' ORA: ',F10.5,' LAT:',F10.5,' LONG:',F10.5 00001200
+ 1 ,' PARAMETRI INTERNI:',2E12.6,I5) 00001210
+ WRITE(11,4000) 00001220
+ 4000 FORMAT(/,45X,' DETTAGLI DEL MOTO DEI PIANETI',/) 00001230
+ WRITE(2,5000) 00001240
+ 5000 FORMAT(/,45X,'ASCENSIONE RETTA IN ORE E MINUTI',/) 00001250
+ WRITE(3,6000) 00001260
+ 6000 FORMAT(/,45X,'DECLINAZIONE IN GRADI E PRIMI',/) 00001270
+ WRITE(4,7000) 00001280
+ 7000 FORMAT(/,45X,'COORDINATE GEOCENTRICHE ECLITTICHE', 00001290
+ 1 ' LAMBDA,BETA',/) 00001300
+ WRITE(2,5500) 00001310
+ WRITE(3,5500) 00001320
+ WRITE(4,5500) 00001330
+ 5500 FORMAT(5X,'DATA:',8X,'SOLE',4X,'MERCURIO',4X,'VENERE',7X, 00001340
+ 2 'MARTE',6X,'GIOVE',4X,'SATURNO',5X,'URANO',5X, 00001350
+ 3 'NETTUNO',4X,'PLUTONE',4X,'LUNA',/) 00001360
+ 102 CONTINUE 00001370
+C ---------------------------------------------------------- 00001380
+ TYPE*,' per questa volta ci credo, ora provo a fare i conti.' 00001390
+C Le metto in radianti,come tutti gli angoli nel programma. 00001400
+ ALATR=(INT(ALAT)+(ALAT-INT(ALAT))/60.D0*100.D0)/360.D0*2*PIGREC 00001410
+ ALONGR=(INT(ALONG)+(ALONG-INT(ALONG))/60.D0*100)/24.D0*2*PIGREC 00001420
+C 00001430
+C --------------------------------------------------------------- 00001440
+C Apertura del file dei dati 00001450
+ OPEN(UNIT=99,FILE='CIBO',STATUS='OLD',
+ 1 ORGANIZATION='INDEXED',ACCESS='KEYED',RECORDTYPE='VARIABLE', 00001470
+ 2 FORM='UNFORMATTED',ERR=105,RECL=250,KEY=(1:4:INTEGER)) 00001480
+C Calcolo valore della chiave 00001490
+ KIAV=-(GIORNO+MESE*100+ANNO*10000) 00001500
+ READ(UNIT=99,KEYGE=KIAV,KEYID=0,ERR=105) 00001510
+ 1 KV,TEMPI,((E(J,JJ),J=1,6),JJ=1,9) 00001520
+ CLOSE(UNIT=99,DISP='KEEP') 00001530
+ CALL JOD(TEMPI,GIOIN,MESIN,IANIN) 00001540
+ IF(STMP.LE.0) GO TO 106 00001550
+ WRITE(11,7700)TEMPI,GIOIN,MESIN,IANIN, 00001560
+ 1 (JJ,(E(J,JJ),J=1,6),JJ=1,9) 00001570
+ 7700 FORMAT(' DATI INIZIALI:JD=',E20.15,' DATA:',F6.3,2I5, 00001580
+ 1 ' ELEMENTI ORBITALI: (TERRA PRIMO POSTO)',/ 00001590
+ 2 ' PIANETA',7X,'I',13X,'OMP',13X,'OMG',14X,'ECC',13X,'ANOM', 00001600
+ 3 12X,'A',/,9(/,1X,I5,6E16.7),/) 00001610
+ GOTO106 00001620
+ 105 TYPE*,' AARGH! HO FAME,HO FAME! NON TROVO IL CIBO.DAT! ARGH'00001630
+ TYPE*,' AARGH! DOVE E` IL MIO CIBO? DAMMELO !' 00001640
+ STOP 00001650
+ 106 CONTINUE 00001660
+ DO 20 II=1,9 00001670
+ E(1,II)=GRARAD(E(1,II)) 00001680
+ E(2,II)=GRARAD(E(2,II)) 00001690
+ E(3,II)=GRARAD(E(3,II)) 00001700
+ E(5,II)=GRARAD(E(5,II)) 00001710
+ 20 CONTINUE 00001720
+C ----------------------------------------------------------------00001730
+C Giorni giuliani interi dal 1900 per il calcolo del tempo sideral00001740
+C e dell'obliquita` eclittica (Meeu00001750
+ T1=DJ(ANNO,MESE,GIORNO) 00001760
+ T1=(T1-2415020.D0)/36525.D0 00001770
+C Obliquita' dell'eclittica dal Meeus 00001780
+ EPS=23.452294-0.0130125*T1-0.00000164*T1**2 00001790
+ EPS=GRARAD(EPS) 00001800
+C Longitudine in ore 00001810
+ ALONGO=ALONGR/(2.D0*PIGREC)*24.D0 00001820
+ ORAD=INT(ORA)+(ORA-INT(ORA))/60.D0*100.D0 00001830
+C ----------------------------------------------------- 00001840
+C tempo siderale finale locale 00001850
+C formule dal libro di Meeus gia citato 00001860
+ TS=6.646056D0+2400.051262*T1+0.00002581*T1**2 00001870
+ TS=MOD(TS,24.D0) 00001880
+ TS=TS+ALONGO+(ORAD)*1.0027379D0 00001890
+ TS=MOD(TS,24.D0) 00001900
+ IF(STMP.GT.0)WRITE(11,*)' tempo siderale:',TS 00001910
+C tempo siderale in radianti 00001920
+ TS=TS/24*2*PIGREC 00001930
+C ----------------------------------------------------------------00001940
+ TFIN=DJ(ANNO,MESE,GIORNO) 00001950
+ TFIN=TFIN+GIOR(ORA) 00001960
+ TINIZ=TEMPI 00001970
+C Metto i tempi in secondi giuliani,per coerenza con G, in S,Kg,Km00001980
+ TFIN=TFIN*86400 00001990
+ TINIZ=TINIZ*86400 00002000
+ TEMPO=TFIN-TINIZ 00002010
+C --------------------------------------------------------------- 00002020
+C intervallo per il calcolo dei pianeti: 36525 giorni.(Fa 1 solo s00002030
+ DTP=86400*36525.D0 00002040
+C correggo l'intervallo in modo che alla fine del loop si 00002050
+C arrivi esattamete al tempo finale 00002060
+ NT=ABS(TEMPO/DTP) 00002070
+ IF(NT.EQ.0) NT=1 00002080
+ DTP=TEMPO/NT ! INSERISCE IL RESTO DELLA DIVISIONE 00002090
+C 00002100
+C Input esplicito del passo 00002110
+ IF(PASSOP.NE.0) DTP=PASSOP*86400 00002120
+C 00002130
+ IF(STMP.GT.0) WRITE(11,7800)TINIZ,TFIN,TEMPO,DTP 00002140
+ 7800 FORMAT(' TEMPO INIZIALE:',F15.2,' TEMPO FINALE:',F15.2, 00002150
+ 1 ' INTERVALLO:',F15.2,' DTP:',E15.9) 00002160
+C ----------------------------------------------------------------00002170
+C ----------------------------------------------------------------00002180
+ T=TINIZ 00002190
+ T=T-DTP 00002200
+C 00002210
+C --------------------------------LOOP SU ORBITE IMPERTURBATE-- 00002220
+C ----------------------------------------------------------------00002230
+ 51 CONTINUE 00002240
+ T=T+DTP 00002250
+C 00002260
+C evoluzione della luna fino al tempo T 00002270
+ CALL LUNAS(XL,R,T/86400.D0,ALAMBL,BETL,PIL) 00002280
+C evoluzione imperturbata dei pianeti su ellissi fino al tempo T 00002290
+ CALL TEMA(X,F,AE,AME,E,T,TINIZ) 00002300
+C Coordinate eliocentriche della luna: 00002310
+ DO 30 J=1,3 00002320
+ 30 X(J,10)=XL(J)+X(J,1) 00002330
+C eventuale rettificazione delle orbite 00002340
+ IF(KRETT.EQ.0) GOTO49 00002350
+ IF(TINIZ.EQ.0) GOTO 49 ! PER IL TEMPO INIZIALE NON RETTIFICA OVVI00002360
+ CALL RETT(E,X,F,DTP,AE) 00002370
+ TINIZ=T ! OPPURE SI USA N FINALE SEMPRE INVECE CAMBIA E' N(T) 00002380
+ DO 48 J=1,9 00002390
+ 48 E(5,J)=AME(J)+E(5,J) 00002400
+C 00002410
+ 49 CONTINUE 00002420
+C ----------------------------------------------------------- 00002430
+C MANCA CORREZIONE AL BARICENTRO TERRA-LUNA !!!!!!!!!!!!!!! 00002440
+C PARE PICCOLA. 00002450
+C ----------------------------------------------------------- 00002460
+ IF(STMP.GT.0)CALL CORDI(T,X,X1,XL,BETA,LAMBDA,ALFA,DELTA,ALAMBL, 00002470
+ 1 BETL) 00002480
+ 50 CONTINUE 00002490
+ IF(T.GE.TFIN.AND.DTP.GE.0.D0.OR.T.LE.TFIN.AND.DTP.LE.0.D0)goto50000002500
+ GO TO 51 00002510
+ 500 CONTINUE 00002520
+C ----------------------------------------------------------------00002530
+c ********* fine loop sugli intervalli di tempo *** 00002540
+C --------------------------------------------------------------- 00002550
+ IF(PASSOP.NE.0) TFIN=T 00002560
+ CALL CORDI(TFIN,X,X1,XL,BETA,LAMBDA,ALFA,DELTA,ALAMBL,BETL) 00002570
+C pianeti e luna nei segni 00002580
+ DO 80 I=1,10 00002590
+ SEGNI(I)=INT(LAMBDA(I)/(2*PIGREC)*12+1) 00002600
+ IF(SEGNI(I).EQ.13) SEGNI(I)=12 00002610
+ 80 CONTINUE 00002620
+C ------------------------------------------------------- 00002630
+C riporto il tempo da secondi a giorni giuliani 00002640
+ TFIN=TFIN/86400. 00002650
+C ------------------------------------------------------- 00002660
+C Calcolo latitudine relativa alle cuspidi delle case 00002670
+ CALL DOMI(CUSP,TS,ALATR,ALONGR) 00002680
+C ------------------------------------------------------- 00002690
+C Metto i pianeti nelle case 00002700
+C 00002710
+ DO 90 I=1,10 00002720
+ DO 95 J=1,12 00002730
+ J1=J+1 00002740
+ J0=J-1 00002750
+ IF(J1.GT.12)J1=1 00002760
+ IF(J0.LT.1) J0=12 00002770
+ IF(LAMBDA(I).GT.CUSP(J).AND.LAMBDA(I).LT.CUSP(J1))GOTO180 00002780
+ IF(CUSP(J).LT.CUSP(J0).AND.LAMBDA(I).LT.CUSP(J))GOTO180 00002790
+ 95 CONTINUE 00002800
+ CASE(I)=12 00002810
+ GO TO 90 00002820
+ 180 CASE(I)=J 00002830
+ 90 CONTINUE 00002840
+ DO 96 I=1,12 00002850
+ ISEGNC(I)=CUSP(I)*12/(2.D0*PIGREC)+1 00002860
+ IF(ISEGNC(I).GT.12) ISEGNC(I)=12 00002870
+ 96 CONTINUE 00002880
+C ------------------------------------------------------- 00002890
+C Stampa del tema natale 00002900
+ CALL STAMPAT(SEGNI,LAMBDA,CUSP,CASE,NOME,ALAT,ALONG,TFIN, 00002910
+ 1 ISEGNC) 00002920
+ STOP 00002930
+ END 00002940
+ SUBROUTINE CORDI(T,X,X1,XL,BETA,LAMBDA,ALFA,DELTA,ALAMBL,BETL) 00002950
+C ************************************************************* 00002960
+C CALCOLA COORDINATE DI TUTTI I TIPI E FA UNA MAREA DI STAMPE 00002970
+C ************************************************************ 00002980
+ IMPLICIT REAL*8 (A-H,O-Z) 00002990
+ COMMON/MASS/MS,MT,ML,AMS(9),G,PIGREC,UA,EPS,PASSOP, 00003000
+ 1 STMP,KRETT 00003010
+ DIMENSION X(3,10),X1(3,10),XL(3),BETA(10),LAMBDA(10) 00003020
+ DIMENSION ALFA(10),DELTA(10),BET1(10),ALAM1(10) 00003030
+ REAL*8 LAMBDA,MS,MT,ML 00003040
+ DIMENSION IGR(10),APR(10),IOR(10),AMINU(10) 00003050
+C 00003060
+ T1=T/86400 00003070
+ CALL JOD(T1,GIORNO,MESE,IANNO) 00003080
+C 00003090
+ IF(STMP.GT.120)WRITE(11,1000)((X(J1,J2),J1=1,3) 00003100
+ 1 ,IANNO,MESE,GIORNO ,J2=1,10) 00003110
+ 1000 FORMAT(' COORDINATE ELIOCENTRICHE RETTANGOLARI PIANETI', 00003120
+ 1 /,(3(1X,E16.10),' DATA:',I6,I5,F6.2)) 00003130
+C -------------------------------------------------------- 00003140
+C calcolo coordinate geocentriche rettangolari dei pianeti 00003150
+ DO 60 J=1,3 00003160
+ 60 X1(J,1)=-X(J,1) 00003170
+ DO 70 I=2,9 00003180
+ DO 70 J=1,3 00003190
+ 70 X1(J,I)=X(J,I)-X(J,1) 00003200
+C 00003210
+C ------------------------------------------------------- 00003220
+C mette neL vettore coordinate la luna all'ultimo posto 00003230
+ X1(1,10)=XL(1) 00003240
+ X1(2,10)=XL(2) 00003250
+ X1(3,10)=XL(3) 00003260
+C 00003270
+ IF(STMP.GT.125)WRITE(11, 2000)((X1(J1,J2),J1=1,3) 00003280
+ 1 ,IANNO,MESE,GIORNO,J2=1,10) 00003290
+ 2000 FORMAT(' COORDINATE GEOCENTRICHE RETTANGOLARI PIANETI', 00003300
+ 1 /,(3(1X,E16.10),' DATA:',I6,I5,F6.2)) 00003310
+C 00003320
+C --------------------------------------------------------- 00003330
+C longitudine eclittica e declinazione pianeti e luna 00003340
+ CALL RETTECL(X1,BETA,LAMBDA) 00003350
+C ------------------------------------------------------ 00003360
+ DO 72 J=1,10 00003370
+ LAMBDA(J)=MOD(LAMBDA(J),2.D0*PIGREC) 00003380
+ BETA(J)=MOD(BETA(J),2.D0*PIGREC) 00003390
+ BET1(J)=RADGRAD(BETA(J)) 00003400
+ ALAM1(J)=RADGRAD(LAMBDA(J)) 00003410
+ 72 CONTINUE 00003420
+ BET1(10)=RADGRAD(BETL) ! DOVREBBERO AD OGNI MODO ESSERE EGUALI 00003430
+ ALAM1(10)=RADGRAD(ALAMBL) 00003440
+ IF(STMP.GT.0) 00003450
+ 1 WRITE(4,1100) GIORNO,MESE,IANNO,(ALAM1(J),BET1(J),J=1,10) 00003460
+ 1100 FORMAT(1X,F5.2,I3,I5,10(1H*,2F5.1)) 00003470
+C ------------------------------------------------------- 00003480
+C ascensione retta e declinazione 00003490
+ CALL ECLEQ(EPS,BETA,LAMBDA,ALFA,DELTA) 00003500
+C --------------------------------------------------------- 00003510
+ DO 75 J=1,10 00003520
+C ALfa in ore e minuti 00003530
+ AAA=ALFA(J)/(2.D0*PIGREC)*24.D0 00003540
+ AMINU(J)=(AAA-INT(AAA))*60.D0 00003550
+ IOR(J)=INT(AAA) 00003560
+C Delta in gradi e primi 00003570
+ AAA=RADGRAD(DELTA(J)) 00003580
+ IGR(J)=INT(AAA) 00003590
+ APR(J)=(AAA-IGR(J))*60.D0 00003600
+ 75 CONTINUE 00003610
+C 00003620
+ IF(STMP.GT.0) WRITE(2,4000)GIORNO,MESE,IANNO, 00003630
+ 1 (IOR(J),AMINU(J),J=1,10) 00003640
+ 4000 FORMAT(1X,F5.2,I3,I5,(10(1H*,I3,1X,F6.2))) 00003650
+ IF(STMP.GT.0) WRITE(3,5000)GIORNO,MESE,IANNO, 00003660
+ 1 (IGR(J),APR(J),J=1,10) 00003670
+ 5000 FORMAT(1X,F5.2,I3,I5,10(1H*,I4,F6.2)) 00003680
+ RETURN 00003690
+ END 00003700
+ SUBROUTINE ANOM(E,MP,MS,G,T,AME,AE,F,T0,PIGREC) 00003710
+C *************************************************** 00003720
+C CALCOLO ANOMALIA DEL PIANETA AL TEMPO T 00003730
+C *************************************************** 00003740
+ IMPLICIT REAL*8 (A-H,O-Z) 00003750
+ REAL*8MS,MP 00003760
+ DIMENSION E(6) 00003770
+ ENNE=SQRT(G*(MP+MS))/SQRT(E(6))**3 00003780
+ AME=ENNE*(T-T0)+E(5) 00003790
+ AE=AME 00003800
+ N=0 00003810
+ AMEG=RADGRAD(AE) 00003820
+ 90 AE1=AME+E(4)*SIN(AE) 00003830
+ IF(ABS(AE1-AE).LT.1.E-10) GO TO80 00003840
+ AE=AE1 00003850
+ N=N+1 00003860
+ IF(N.GT.100) GOTO 80 00003870
+ GO TO 90 00003880
+ 80 CONTINUE 00003890
+ AMEG1=RADGRAD(AE) 00003900
+ W=SQRT((1+E(4))/(1-E(4)))*TAN(AE/2) 00003910
+ F=ATAN(W)*2 00003920
+C marchingegno per cacciare nel giusto quadrante la 00003930
+C arcotangente che mette fra 90 e - 90 00003940
+ FXP=ABS(F-AE) 00003950
+ IF(FXP.GT.PIGREC/2.AND.FXP.LT.PIGREC*3.D0/2)F=F+PIGREC 00003960
+ F=MOD(F,2.D0*PIGREC) 00003970
+ RETURN 00003980
+ END 00003990
+ SUBROUTINE TEMA(X,F,AE,AME,E,T,T0) 00004000
+C **************************************************** 00004010
+C EVOLUZIONE DEI PIANETI SU ORBITE IMPERTURBATE 00004020
+C **************************************************** 00004030
+ IMPLICIT REAL*8 (A-H,O-Z) 00004040
+ REAL*8 MS,MT,ML 00004050
+ COMMON/MASS/MS,MT,ML,AM(9),G,PIGREC,UA,EPS,PASSOP, 00004060
+ 1 STMP,KRETT 00004070
+ DIMENSION X(3,9),F(9),AE(9),E(6,9) 00004080
+ DIMENSION ANM(3,9),AME(9) 00004090
+ CALL JOD(T/86400.D0,GIORNO,MESE,IANNO) 00004100
+C 00004110
+ TYPE*,' MUOVO I PIANETI !',' DATA:',GIORNO,MESE,IANNO 00004120
+C 00004130
+ DO 10 I=1,9 00004140
+C IF(STMP.GT.100.D0)WRITE(11,1000)I,T,T0,GIORNO,MESE,IANNO 00004150
+C1000 FORMAT(' PIANETA',I5,' TFINALE:',E20.10,' TINIZIAL:',E20.10,00004160
+C 1 ' DATA: ',F6.2,2I5) 00004170
+C calcolo anomalia vera e raggio vettore 00004180
+ CALL ANOM(E(1,I),AM(I),MS,G,T,AME(I),AE(I),F(I),T0,PIGREC) 00004190
+ R=E(6,I)*(1-E(4,I)*COS(AE(I))) 00004200
+C calcolo coordinate eliocentriche rettangolari eclittiche 00004210
+ CALL NODECL(E(1,I),X(1,I),R*COS(E(2,I)+F(I)),R*SIN(E(2,I)+F(I)) 00004220
+ 1 ,0.D0) 00004230
+ 10 CONTINUE 00004240
+ IF(STMP.LE.50.D0) GOTO100 00004250
+ DO 20 I=1,9 00004260
+ ANM(1,I)=RADGRAD(AME(I)) 00004270
+ ANM(2,I)=RADGRAD(AE(I)) 00004280
+ ANM(3,I)=RADGRAD(F(I)) 00004290
+ 20 CONTINUE 00004300
+ IF(STMP.GT.0) WRITE(11,3000) 00004310
+ 1 ((IANNO,MESE,GIORNO,JJ,(ANM(J,JJ),J=1,3)),JJ=1,9) 00004320
+ 3000 FORMAT(' DATA:',I5,I5,F6.2,' PIANETA:',I5,' ANOMALIA MEDIA:'00004330
+ 1 ,F8.3,' ANOMALIA ECCENTRICA:',F8.3,' VERA:',F8.3) 00004340
+ 100 CONTINUE 00004350
+ RETURN 00004360
+ END 00004370
+ SUBROUTINE LUNAS(XL,R,TT,ALAMBD,BET,PI) 00004380
+C **************************************************** 00004390
+C MOTO DELLA LUNA CON LE FORMULE PRESE DAL LIBRO DI 00004400
+C MEEUS - ASTRONOMICAL FORMULAE FOR CALCULATORS 00004410
+C **************************************************** 00004420
+ IMPLICIT REAL*8 (A-H,O-Z) 00004430
+ REAL*8 M,M1,L1 00004440
+ DIMENSION XL(3) 00004450
+ T=(TT-2415020D0)/36525.D0 00004460
+C Tempo in secoli giuliani da 0.5 GENN 1900 00004470
+C TERMINI SVILUPPI PERTURBATIVI 00004480
+ T2=T*T 00004490
+ L1=270.434164+481267.8831*T-0.001133*T2 00004500
+ M=358.475833+35999.0498*T-0.00015*T2 00004510
+ M1=296.104608+477198.8491*T+0.009192*T2 00004520
+ D=350.737486+445267.1142*T-0.001436*T2 00004530
+ F=11.250889+483202.0251*T-0.003211*T2 00004540
+ E=1-0.002495*T-0.00000752*T2 00004550
+ E2=E*E 00004560
+ L1=MOD(L1,360.D0) 00004570
+ M=MOD(M,360.D0) 00004580
+ M1=MOD(M1,360.D0) 00004590
+ D=MOD(D,360.D0) 00004600
+ F=MOD(F,360.D0) 00004610
+C 00004620
+C 00004630
+C TYPE*,'L1,M.M1,D,F,T,E,T2',L1,M,M1,D,F,T,E,T2 00004640
+ M=GRARAD(M) 00004650
+ M1=GRARAD(M1) 00004660
+ D=GRARAD(D) 00004670
+ F=GRARAD(F) 00004680
+C Longitudine eclittica 00004690
+ AL=E2*(0.002249*SIN(2*D-2*M)-0.002079*SIN(2*M)+ 00004700
+ 1 0.002059*SIN(2*D-M1-2*M)+0.000717*SIN(M1-2*M) + 00004710
+ 2 0.000704*SIN(M1-2*M-2*D) ) 00004720
+ AL=AL+E*(-0.185596*SIN(M)+0.057212*SIN(2*D-M-M1)+ 00004730
+ 2 0.045874*SIN(2*D-M)+0.041024*SIN(M1-M)- 00004740
+ 3 0.030465*SIN(M+M1)-0.007910*SIN(M-M1+2*D)- 00004750
+ 4 0.006783*SIN(2*D+M)+0.005*SIN(M+D)+ 00004760
+ 5 0.004049*SIN(M1-M+2*D)+0.002695*SIN(2*M1-M)+ 00004770
+ 6 0.002396*SIN(2*D-M-2*M1)-0.002125*SIN(2*M1+M)+ 00004780
+ 7 0.00122*SIN(4*D-M-M1)-0.000811*SIN(M+M1+2*D)+ 00004790
+ 8 0.000761*SIN(4*D-M-2*M1)+0.000693*SIN(M-2*M1+2*D)+ 00004800
+ 9 0.000598*SIN(2*D-M-2*F)+0.000521*SIN(4*D-M) ) 00004810
+ AL=AL+L1+6.28875*SIN(M1)+1.274018*SIN(2*D-M1)+0.658309*SIN(2*d)+ 00004820
+ 1 0.213616*SIN(2*M1)-0.114336*SIN(2*F)+ 00004830
+ 2 0.058793*SIN(2*D-2*M1)+0.053320*SIN(2*D+M1)- 00004840
+ 3 0.034718*SIN(D)+0.015326*SIN(2*D-2*F)-0.012528*SIN(2*F+M1)-00004850
+ 4 0.01098*SIN(2*F-M1)+0.010674*SIN(4*D-M1)+0.010034*SIN(3*M1)+00004860
+ 5 0.008548*SIN(4*D-2*M1)+0.005162*SIN(M1-D)+ 00004870
+ 6 0.003996*SIN(2*M1+2*D)+0.003862*SIN(4*D)+ 00004880
+ 7 0.003665*SIN(2*D-3*M1)+0.002602*SIN(M1-2*F-2*D)- 00004890
+ 8 0.002349*SIN(M1+D)-0.001773*SIN(M1+2*D-2*F)- 00004900
+ 9 0.001595*SIN(2*F+2*D)-0.00111*SIN(2*M1+2*F)+ 00004910
+ A 0.000892*SIN(M1-3*D)+0.00055*SIN(M1+4*D)+ 00004920
+ B 0.000538*SIN(4*M1)+0.000486*SIN(2*M1-D) 00004930
+C Latitidine eclitt00004940
+ B=E2*0.000306*SIN(2*D-2*M-F) 00004950
+ B=B+E*(0.008247*SIN(2*D-M-F)+0.003372*SIN(F-M-2*D)+ 00004960
+ 1 0.002472*SIN(2*D+F-M-M1)+0.002222*SIN(2*D+F-M)+ 00004970
+ 2 0.002072*SIN(2*D-F-M-M1)+0.001877*SIN(F-M+M1)- 00004980
+ 3 0.001803*SIN(F+M)+ 00004990
+ 4 0.00157*SIN(M1-M-F)-0.001481*SIN(F+M+M1)+ 00005000
+ 5 0.001417*SIN(F-M-M1)+0.00135*SIN(F-M)+ 00005010
+ 6 0.000492*SIN(2*D+M1-M-F)-0.000367*SIN(M+F+2*D-M1)- 00005020
+ 7 0.000353*SIN(M+F+2*D)+0.000317*SIN(2*D+F-M+M1) ) 00005030
+ B=B+5.128189*SIN(F)+0.280606*SIN(M1+F)+0.277693*SIN(M1-F)+ 00005040
+ 1 0.173238*SIN(2*D-F)+0.055413*SIN(2*D+F-M1)+ 00005050
+ 2 0.046272*SIN(2*D-F-M1)+0.032573*SIN(2*D+F)+ 00005060
+ 3 0.017198*SIN(2*M1+F)+0.009267*SIN(2*D+M1-F)+ 00005070
+ 4 0.008823*SIN(2*M1-F)+0.004323*SIN(2*D-F-2*M1)+ 00005080
+ 5 0.0042*SIN(2*D+F+M1)+0.001828*SIN(4*D-F-M1)- 00005090
+ 6 0.00175*SIN(3*F)-0.001487*SIN(F+D)+0.00133*SIN(F-D)+ 00005100
+ 7 0.001106*SIN(F+3*M1)+0.00102*SIN(4*D-F)+ 00005110
+ 8 0.000833*SIN(F+4*D-M1)+0.000781*SIN(M1-3*F)+ 00005120
+ 9 0.00067*SIN(F+4*D-2*M1)+0.000606*SIN(2*D-3*F)+ 00005130
+ A 0.000597*SIN(2*D+2*M1-F)+0.00045*SIN(2*M1-F-2*D)+ 00005140
+ B 0.000439*SIN(3*M1-F)+0.000423*SIN(F+2*D+2*M1)+ 00005150
+ C 0.000422*SIN(2*D-F-3*M1)+0.000331*SIN(F+4*D)- 00005160
+ D 0.000283*SIN(M1+3*F) 00005170
+C 00005180
+ PI=E2*0.000026*COS(2*D-2*M) 00005190
+ PI=PI+E*(0.000533*COS(2*D-M)+0.000401*COS(2*D-M-M1)+ 00005200
+ 1 0.000320*COS(M1-M)-0.000264*COS(M+M1)-0.000111*COS(M)- 00005210
+ 2 0.000083*COS(2*D+M)+0.000064*COS(2*D-M+M1)- 00005220
+ 2 0.000063*COS(2*D+M-M1)+ 00005230
+ 3 0.000041*COS(M+D)+0.000035*COS(2*M1-M)- 00005240
+ 4 0.000029*COS(2*M1+M)+0.000019*COS(4*D-M-M1) ) 00005250
+ PI=PI+0.9507240+0.051818*COS(M1)+0.009531*COS(2*D-M1)+ 00005260
+ 1 0.007843*COS(2*D)+0.002824*COS(2*M1)+ 00005270
+ 2 0.000857*COS(2*D+M1)-0.000271*COS(D)- 00005280
+ 3 0.000198*COS(2*F-M1)+0.000173*COS(3*M1)+ 00005290
+ 4 0.000167*COS(4*D-M1)+0.000103*COS(4*D-2*M1)- 00005300
+ 5 0.000084*COS(2*M1-2*D)+0.000079*COS(2*D+2*M1)+ 00005310
+ 6 0.000072*COS(4*D)-0.000033*COS(3*M1-2*D)- 00005320
+ 7 0.00003*COS(M1+D)-0.000029*COS(2*F-2*D)- 00005330
+ 8 0.000023*COS(2*F-2*D+M1) 00005340
+C 00005350
+ PI=GRARAD(PI) 00005360
+ BET=GRARAD(B) 00005370
+ ALAMBD=GRARAD(AL) 00005380
+C 00005390
+C Raggio vettore 00005400
+ R=6378.14/SIN(PI) 00005410
+C 00005420
+C Coordinate Rettangolari Geocentriche Eclittiche 00005430
+ XL(1)=R*COS(BET)*COS(ALAMBD) 00005440
+ XL(2)=R*COS(BET)*SIN(ALAMBD) 00005450
+ XL(3)=R*SIN(BET) 00005460
+C 00005470
+ RETURN 00005480
+ END 00005490
+ SUBROUTINE RETT(E,X,F,DT,AE) 00005500
+C *********************************************************** 00005510
+C FA UNO STEP DELLA SOLUZIONE COL METODO DI TAYLOR DELLE 00005520
+C EQUAZIONI DI LAGRANGE NELLA FORMA DI GAUSS 00005530
+C ********************************************************* 00005540
+ IMPLICIT REAL*8 (A-H,O-Z) 00005550
+ COMMON/MASS/AMS,AMT,AML,AM(9),G,PIGREC,UA,EPS,PASSOP, 00005560
+ 1 STMP,KRETT 00005570
+ DIMENSION F(9),S(3),FOR(3),E(6,9),X(3,9),AE(9) 00005580
+C E=elementi orbitali 00005590
+C X=posizioni dei pianeti 00005600
+C AMS=massa del sole 00005610
+C AM=masse pianeti 00005620
+C F=anomalie vere dei pianeti 00005630
+C G=costante gravitazionale 00005640
+C FOR=forza perturbante nel sistema eclittico 00005650
+C S=coordinate della forza perturbante nel piano orbitale 00005660
+C AE=anomalia eccentrica 00005670
+C --------------------------------------------------------------- 00005680
+ TYPE*,' RETTIFICO UN POCO DI ORBITE,NON SI SA MAI', 00005690
+ 1 ' CHE TORNI BUONO' 00005700
+C loop sui pianeti ---------------------------------------------- 00005710
+ DO 10 I=1,9 00005720
+C calcolo forza perturbante 00005730
+ FOR(1)=0. 00005740
+ FOR(2)=0. 00005750
+ FOR(3)=0. 00005760
+ DO 20 II=1,9 00005770
+C salto il pianeta stesso,che non si autoperturba 00005780
+ IF(II.EQ.I) GOTO 20 00005790
+ RI3=(SQRT(X(1,II)**2+X(2,II)**2+X(3,II)**2))**3 00005800
+ RRI3=(SQRT((X(1,I)-X(1,II))**2+(X(2,I)-X(2,II))**2+ 00005810
+ 1 (X(3,I)-X(3,II))**2))**3 00005820
+ DO 30 IJ=1,3 00005830
+ 30 FOR(IJ)=FOR(IJ)-G*AM(II)*((X(IJ,I)-X(IJ,II))/RRI3+X(IJ,II)/RI3) 00005840
+C IF(STMP.GT.150) WRITE(11,*)'FORZA PERTURBANTE SUL PIANETA:00005850
+C 1 FOR,'COMPRESO IL PIANETA:',II 00005860
+ 20 CONTINUE 00005870
+C ------------------------------------------------------------- 00005880
+C conversione a coordinate sul piano orbitale 00005890
+ CALL ECLNOD(E(1,I),F(I),S,FOR) 00005900
+ IF(STMP.GT.150.D0) WRITE(11,1000) S,I,DT 00005910
+ 1000 FORMAT(' FORZA PERTURBANTE TOTALE SU PIANO ORBITA ',3E15.7, 00005920
+ 1 ' PIANETA:',I3,' TEMPO:',E12.6) 00005930
+C ------------------------------------------------------------ 00005940
+C calcolo pezzi formule successive 00005950
+ AMI=G*(AMS+AM(I)) 00005960
+ ENNE=SQRT(AMI)/SQRT(E(6,I))**3 00005970
+ SINF=SIN(F(I)) 00005980
+ COSF=COS(F(I)) 00005990
+ R=SQRT(X(1,I)**2+X(2,I)**2+X(3,I)**2) 00006000
+ P=E(6,I)*(1-E(4,I)**2) 00006010
+ COSU=COS(E(2,I)+F(I)) 00006020
+ SINU=SIN(E(2,I)+F(I)) 00006030
+ UNME=SQRT(1-E(4,I)**2) 00006040
+ SIN2I2=SIN(E(1,I)/2)**2 00006050
+ A2=E(6,I)**2 00006060
+C --------------------------------------------------------- 00006070
+C calcolo DA 00006080
+ DADT=2./(ENNE*UNME)*(S(1)*E(4,I)*SINF+P*S(2)/R) 00006090
+ DA=DADT*DT 00006100
+C calcolo DE 00006110
+ DE=(UNME/(ENNE*E(6,I))*(S(1)*SINF+S(2)*(COS(AE(I))+COSF)))*DT 00006120
+C calcolo DI 00006130
+ DI=(S(3)*R*COSU/(ENNE*A2*UNME))*DT 00006140
+C calcolo OG 00006150
+ IF(E(1,I).NE.0.D0) THEN 00006160
+ DOGDT=(S(3)*R*SINU/(ENNE*A2*UNME*SIN(E(1,I)))) 00006170
+ ELSE 00006180
+ DOGDT=0.D0 00006190
+ END IF 00006200
+ DOG=DOGDT*DT 00006210
+C calcolo DOS 00006220
+ DOSDT=UNME/(ENNE*E(6,I)*E(4,I))* 00006230
+ 1 (-S(1)*COSF+S(2)*(1+R/P)*SINF)+2*DOGDT*SIN2I2 00006240
+ DOS=DOSDT*DT 00006250
+C calcolo DEP 00006260
+ DEPDT=DOSDT*E(4,I)**2/(1+UNME)+2*DOGDT*UNME*SIN2I2-2*R*S(1)/ 00006270
+ 1 (ENNE*A2) 00006280
+ DEP=DEPDT*DT 00006290
+C ------------------------------------------------------------ 00006300
+C calcolo DOP 00006310
+ DOP=DOS-DOG 00006320
+C ----------------------------------------------------------- 00006330
+C correggo gli elementi dell'orbita 00006340
+ E(1,I)=E(1,I)+DI 00006350
+ E(2,I)=E(2,I)+DOP 00006360
+ E(3,I)=E(3,I)+DOG 00006370
+ E(4,I)=E(4,I)+DE 00006380
+ E(5,I)=DEP 00006390
+ E(6,I)=E(6,I)+DA 00006400
+ IF(STMP.LE.0.D0) GOTO10 00006410
+ DIS1=RADGRAD(DI) 00006420
+ DOPS1=RADGRAD(DOP) 00006430
+ DOGS1=RADGRAD(DOG) 00006440
+ DELLE1=RADGRAD(DEP) 00006450
+ IF(STMP.GT.0) WRITE(11,6000) 00006460
+ 1 DIS1,DOPS1,DOGS1,DE,DELLE1,DA 00006470
+ 6000 FORMAT(' CORREZIONI A:DI,DOP,DOG,DE,DELLE,DA',8E15.5) 00006480
+ 10 CONTINUE 00006490
+ RETURN 00006500
+ END 00006510
+ SUBROUTINE DOMI(CUSPL,T,LAT,LONG) 00006520
+C ******************************************************* 00006530
+C CALCOLA LA LONGITUDINE ECLITTICA DELLE CUSPIDI 00006540
+C DELLE 12 CASE. 00006550
+C Sembra che questo metodo sia quello di un certo monaco 00006560
+C Placido , risalente al 1650 circa. Costui, individuati 00006570
+C ascendente e discendente (intersezioni fra eclittica 00006580
+C ed orizzonte) e medium ed imum coeli ( intersezioni 00006590
+C fra meridiano del luogo ed eclittica ); 00006600
+C divide in tre parti gli archi delimitati 00006610
+C sull'equatore dai meridiani di questi 4 punti 00006620
+C Individua cosi' sull'equatore 12 punti i cui meridiani 00006630
+C tagliano a fette l'eclittica; le fette sono le case. 00006640
+C Oltre questo modo di calcolare le case ne esistono 00006650
+C altri 4 o 5, ovviamente piu' semplici, ma questo va 00006660
+C di moda fra gli astrologi perche' loro mica fanno i 00006670
+C conti; vanno a vedere sulle tavole e un certo Rafael 00006680
+C calcola tavole molto diffuse con questo sistema. 00006690
+C Forse non ci abbiamo preso a riprodurre il metodo di 00006700
+C Placido (cosa estremamente probabile visto che lo 00006710
+C abbiamo dedotto da libri di astrologia e come noto 00006720
+C gli astrologi moderni per lo piu' non sanno contare) 00006730
+C In questo caso abbiamo inventato un metodo nuovo, 00006740
+C la cui verita' assoluta sara' affermata in base ad 00006750
+C argomenti inoppugnabili basati sulla scienza, la 00006760
+C cabala, un po' di numerologia, sedute spiritiche ed 00006770
+C altro,come di prassi. 00006780
+C --------------------------------------------------------- 00006790
+ IMPLICIT REAL*8 (A-H,O-Z) 00006800
+ COMMON/MASS/MS,MT,ML,AMS(9),G,PIGREC,UA,EPS,PASSOP, 00006810
+ 1 STMP,KRETT 00006820
+ DIMENSION CUSPL(12),A(3,3),B(3,3),C(3,3),MC(3),IC(3) 00006830
+ DIMENSION ASC(3),DISC(3),CUSP(3,12) 00006840
+ DIMENSION V(3),VV(3),ANGOL(12) 00006850
+ REAL*8 MC,IC,LAT,LONG,MS,MT,ML 00006860
+ TYPE*,' LE CASE, LE CASE! DOMIFICO!' 00006870
+C trasformazione con cui 00006880
+C metto asse x nel primo verticale (rotazione sul piano equatorial00006890
+ COST=COS(T) 00006900
+ SINT=SIN(T) 00006910
+ B(1,1)=COST 00006920
+ B(1,2)=SINT 00006930
+ B(1,3)=0.D0 00006940
+ B(2,1)=-SINT 00006950
+ B(2,2)=COST 00006960
+ B(2,3)=0.D0 00006970
+ B(3,1)=0.D0 00006980
+ B(3,2)=0.D0 00006990
+ B(3,3)=1.D0 00007000
+C definisco medium coeli. Nel sistema dell'equatore, con l'asse x 00007010
+C lungo il meridiano del luogo, si tratta del vettore 1,0,0 . 00007020
+C Lo ruoto nell'equatore in modo da avere x lungo gamma applicando00007030
+C la traformazione inversa di B. 00007040
+C SI PUO' FARE ANCHE COME TRASFORMAZIONE DI ANGOLO DA COORDINATE 00007050
+C ALTAZIMUTALI IN ECLITTICHE! 00007060
+ MC(1)=1.D0 00007070
+ MC(2)=0.D0 00007080
+ MC(3)=0.D0 00007090
+ CALL trasf(ic,b) 00007100
+C L'imum coeli e' all'oppsto del medium coeli (riflessione del vet00007110
+ IC(1)=-MC(1) 00007120
+ IC(2)=-MC(2) 00007130
+ IC(3)=-MC(3) 00007140
+ TS=T-PIGREC/2. 00007150
+C Definisco l'ascendente con la risoluzione di un pasticciato 00007160
+C TRAINGOLO SFERICO. SPERIAMO DI AVERCI PRESO QUESTA VOLTA! 00007170
+ TG=TAN(LAT) 00007180
+ IF(TG.GT.0.D0) THEN 00007190
+ ALM=ATAN2(SIN(TS), 00007200
+ 1 COS(TS)*COS(EPS)+SIN(EPS)/TG ) 00007210
+ ELSE 00007220
+ ALM=0.D0 00007230
+ ENDIF 00007240
+ ALF=ATAN2(COS(EPS)*SIN(ALM),-COS(ALM)) 00007250
+ ASC(1)=COS(ALF) 00007260
+ ASC(2)=SIN(ALF) 00007270
+ ASC(3)=0.0 00007280
+C WRITE(6,1000)ALF,ASC(1),ASC(2),ASC(3) 00007290
+ 1000 FORMAT(' COEFF, ASC;',4E20.10) 00007300
+C definisco il discendente come riflesso dell'ascendente (a 180gra00007310
+ DISC(1)=-ASC(1) 00007320
+ DISC(2)=-ASC(2) 00007330
+ DISC(3)=-ASC(3) 00007340
+C definisco l'angolo fra ascendente e punto gamma 00007350
+ V(1)=1.D0 00007360
+ V(2)=0.D0 00007370
+ V(3)=0.D0 00007380
+ THET=ANGOLO(V,ASC) 00007390
+ IF(ASC(2).GT.0) THET=2*PIGREC-THET 00007400
+C WRITE(6,*)'THET',RADGRAD(THET) 00007410
+ ANGOL(1)=THET 00007420
+ D1=ANGOLO(IC,ASC)/3.D0 00007430
+C WRITE(6,*)'D1 D1 D1 D1',RADGRAD(D1) 00007440
+ ANGOL(2)=ANGOL(1)+D1 00007450
+ ANGOL(3)=ANGOL(2)+D1 00007460
+ ANGOL(4)=ANGOL(3)+D1 00007470
+ D1=ANGOLO(IC,DISC)/3.D0 00007480
+C WRITE(6,*)'D1 D1 D1 D1',RADGRAD(D1) 00007490
+ ANGOL(5)=ANGOL(4)+D1 00007500
+ ANGOL(6)=ANGOL(5)+D1 00007510
+ ANGOL(7)=ANGOL(6)+D1 00007520
+ D1=ANGOLO(MC,DISC)/3.D0 00007530
+C WRITE(6,*)'D1 D1 D1 D1',RADGRAD(D1) 00007540
+ ANGOL(8)=ANGOL(7)+D1 00007550
+ ANGOL(9)=ANGOL(8)+D1 00007560
+ ANGOL(10)=ANGOL(9)+D1 00007570
+ D1=ANGOLO(MC,ASC)/3.D0 00007580
+C WRITE(6,*)'D1 D1 D1 D1',RADGRAD(D1) 00007590
+ ANGOL(11)=ANGOL(10)+D1 00007600
+ ANGOL(12)=ANGOL(11)+D1 00007610
+C trasformazione degli angoli da equatoriale ad eclittico. 00007620
+ COSE=COS(EPS) 00007630
+ DO 40 I=1,12 00007640
+ CUSPL(I)=ATAN2(COSE*SIN(ANGOL(I)),COS(ANGOL(I))) 00007650
+ F1112=RADGRAD(CUSPL(I)) 00007660
+ IF (CUSPL(I).LT.0) CUSPL(I)=CUSPL(I)+2.D0*PIGREC 00007670
+ F1113=RADGRAD(CUSPL(I)) 00007680
+ 40 CONTINUE 00007690
+ RETURN 00007700
+ END 00007710
+ SUBROUTINE STAMPAT(SEGNI,LONG,CUSPIDI,CASE,NOME,ALAT,ALONG,T, 00007720
+ 1 ISEGNC) 00007730
+C ******************************************************** 00007740
+C STAMPA DEL TEMA NATALE 00007750
+C ******************************************************* 00007760
+ IMPLICIT REAL*8 (A-H,O-Z) 00007770
+ CHARACTER*80 GLUMESS(3) 00007780
+ DIMENSION SEGNI(10),CASE(10),CUSPIDI(12),LONG(10) 00007790
+ DIMENSION ISEGNC(12) ! SEGNI DELLE CUSPIDI 00007800
+ REAL*8 LONG 00007810
+ CHARACTER*30,NOME 00007820
+ CHARACTER*10,ZODIAC(12),PIANET(10) 00007830
+ INTEGER SEGNI,CASE 00007840
+ DATA PIANET/' SOLE ','MERCURIO ','VENERE ', 00007850
+ 1 'MARTE ','GIOVE ','SATURNO ','URANO ', 00007860
+ 2 'NETTUNO ','PLUTONE ','LUNA '/ 00007870
+ DATA ZODIAC/'ARIETE ','TORO ','GEMELLI ', 00007880
+ 1 'CANCRO ','LEONE ','VERGINE ','BILANCIA ', 00007890
+ 2 'SCORPIONE ','SAGITTARIO','CAPRICORNO','ACQUARIO ', 00007900
+ 3 'PESCI '/ 00007910
+ WRITE(1,*) ' '
+ WRITE(1,*) ' Oroscopo offerto dalla: '
+ WRITE(1,*) ' FRANKENSTEIN BUILDING CORPORATION '
+ WRITE(1,*) ' BOLOGNA-ITALY '
+ WRITE(1,*) ' '
+ WRITE(1,*)'?????????????????????????????????????????????' 00007920
+ WRITE(1,*)' TEMA NATALE DI :',NOME 00007930
+ CALL JOD(T,G,M,JA) 00007940
+ WRITE(1,*)' NATO IL :',G,M,JA 00007950
+ WRITE(1,*)' ALLA LATITUDINE:',ALAT 00007960
+ WRITE(1,*)' ALLA LONGITUDINE:',ALONG 00007970
+ WRITE(1,*)'?????????????????????????????????????????????' 00007980
+C stampe delle modalita' di uso 00007990
+ OPEN(UNIT=98,FILE='GLUKIE',STATUS='OLD'
+ 1 ,ORGANIZATION='INDEXED',ACCESS='KEYED', 00008010
+ 2 RECORDTYPE='VARIABLE',FORM='UNFORMATTED', 00008020
+ 3 KEY=(1:4:INTEGER)) 00008030
+ DO 30 KIAV=10201,10215 00008040
+ READ(UNIT=98,KEYEQ=KIAV) KIAI,GLUMESS 00008050
+ WRITE(1,*) GLUMESS 00008060
+ 30 CONTINUE 00008070
+ WRITE(1,*)' PIANETI NEI SEGNI:' 00008080
+ DO 10 I=1,10 00008090
+ WRITE(1,*)' ',PIANET(I),' IN ',ZODIAC(SEGNI(I)), 00008100
+ 1 ' Long:',RADGRAD(LONG(I)) 00008110
+ 10 CONTINUE 00008120
+ WRITE(1,*) ' PIANETI NELLA CASE:' 00008130
+ DO 20 I=1,10 00008140
+ WRITE(1,*)' ',PIANET(I),' IN CASA:',CASE(I) 00008150
+ 20 CONTINUE 00008160
+ WRITE(1,*)' POSIZIONI DELLE CUSPIDI DELLE 12 CASE ' 00008170
+ WRITE(1,*)' La prima cuspide e'' l''Ascendente,'
+ WRITE(1,*)' la quarta l''Imum Coeli, la settima il Discendente'
+ WRITE(1,*)' la decima il Medium Coeli.'
+ WRITE(1,1000) (J,RADGRAD(CUSPIDI(J)), 00008180
+ 1 ZODIAC(ISEGNC(J)),J=1,12) 00008190
+ 1000 FORMAT( ' CUSPIDE DELLA ',I5,' CASA IN:',F10.5,' SEGNO:',A10) 00008200
+ WRITE(1,*)' FINE DEL TEMA NATALE...................... ' 00008210
+ WRITE(1,*)' UN CUPO DESTINO SI PROFILA ALL''ORIZZONTE?' 00008220
+ WRITE(1,*)' STAI PER COLLASSARE A BUCO NERO?' 00008230
+ WRITE(1,*)' PIOVERA'' DOMANI O VEDRAI STELLE ?' 00008240
+ WRITE(1,*)' ESPLOSIONI DI SUPERNOVAE ALL''ORIZZONTE?' 00008250
+ WRITE(1,*)' MILLE INTERROGATIVI CI ASSILLANO.' 00008260
+ WRITE(1,*)' LA RISPOSTA LA TROVERAI NELLA PARTE INTERPRETATIVA.' 00008270
+ WRITE(1,*)'$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$' 00008280
+ WRITE(1,1500) 00008290
+ 1500 FORMAT(////,' PARTE INTERPRETATIVA: (AUGURI!!!!!!)',///,1x, 00008300
+ 1 80(1H$),//) 00008310
+C Stampa del responso 00008320
+C Lettura del sole 00008330
+ DO 40 I=1,3 00008340
+ ICHIAV=I*10000+100+SEGNI(1) 00008350
+ READ(UNIT=98,KEYEQ=ICHIAV,ERR=555) KIAI,GLUMESS 00008360
+ WRITE(1,*) GLUMESS 00008370
+ 40 CONTINUE 00008380
+C Lettura dei pianeti nei segni e nelle case 00008390
+ DO 50 K=1,10 00008400
+ ICHIAV=40000+K*100+SEGNI(K) 00008410
+ READ(UNIT=98,KEYEQ=ICHIAV,ERR=555) KIAI,GLUMESS 00008420
+ WRITE(1,*) GLUMESS 00008430
+ 50 CONTINUE 00008440
+C Pianeti nelle case 00008450
+ DO 52 K=1,10 00008460
+ ICHIAV=50000+K*100+CASE(K) 00008470
+ READ(UNIT=98,KEYEQ=ICHIAV,ERR=555) KIAI,GLUMESS 00008480
+ WRITE(1,*) GLUMESS 00008490
+ 52 CONTINUE 00008500
+C LETTURA CASE NEI SEGNI 00008510
+ DO 60 J=1,12 00008520
+ ICHIAV=60000+J*100+ISEGNC(J) 00008530
+ READ(UNIT=98,KEYEQ=ICHIAV,ERR=555) KIAI,GLUMESS 00008540
+ WRITE(1,*) GLUMESS 00008550
+ 60 CONTINUE 00008560
+ CLOSE(UNIT=98,DISP='KEEP') 00008570
+ RETURN 00008580
+ 555 WRITE(1,5555) 00008590
+C 5555 FORMAT(' CI SPIACE MA IL SUO DESTINO E'' MANCHEVOLE MANCANDO00008600
+C 1 un MESSAGGIO INTERPRETATIVO') 00008610
+ 5555 FORMAT( ' L''oroscopo e'' finito!')
+ GO TO 60 00008620
+ END 00008630
+ FUNCTION ANGOLO(X,Y) 00008640
+C ************************************** 00008650
+C CALCOLA L'ANGOLO FRA I VETTORI X ED Y 00008660
+C ****************************************** 00008670
+ IMPLICIT REAL*8 (A-H,O-Z) 00008680
+ DIMENSION X(3),Y(3) 00008690
+ BSX=SQRT(X(1)**2+X(2)**2+X(3)**2) 00008700
+ BSY=SQRT(Y(1)**2+Y(2)**2+Y(3)**2) 00008710
+ SCAL=X(1)*Y(1)+X(2)*Y(2)+X(3)*Y(3) 00008720
+ ANGOLO=ACOS(SCAL/(BSX*BSY)) 00008730
+ RETURN 00008740
+ END 00008750
+ FUNCTION DJ(Y,M,D) 00008760
+C *******************************************
+C Giorno giuliano
+C *******************************************
+ REAL*8 C,D,DJ,A,B,G 00008770
+ INTEGER*4 Y 00008780
+ C=Y+M*1.E-2+D*1.E-4 00008790
+C CALCOLO DELLA CORREZIONE DI GREGORIO 00008800
+ B=0 00008810
+ IF(C.LE.1582.1015) GOTO 100 00008820
+ A=INT(Y/100.) 00008830
+ B=2-A+INT(A/4.) 00008840
+ 100 CONTINUE 00008850
+ IF(M.LE.2)THEN 00008860
+ MY=Y-1 00008870
+ MM=M+12 00008880
+ ELSE 00008890
+ MY=Y 00008900
+ MM=M 00008910
+ END IF 00008920
+ CONTINUE 00008930
+ G=0. 00008940
+ IF(C.LT.0.) G=.75 00008950
+ DJ=INT(365.25*MY-G)+INT(30.6*(MM+1))+D+1720994.5+B 00008960
+C TYPE*,'GIORNO= ',D,' MESE= ',M,' ANNO= ',Y,' DJ=',DJ 00008970
+ RETURN 00008980
+ END 00008990
+ SUBROUTINE ECLEQ(EPS,BETA,LAMBDA,ALFA,DELTA) 00009000
+C ****************************************************** 00009010
+C DA COORDINATE ECLITTICHE AD ASCENSIONE RETTA E DECLINAZIONE 00009020
+C ************************************************************ 00009030
+ IMPLICIT REAL*8 (A-H,O-Z) 00009040
+ DATA PIGREC/3.14159265358979323846/ 00009050
+ REAL*8 BETA(10),LAMBDA(10),ALFA(10),DELTA(10) 00009060
+ COSE=COS(EPS) 00009070
+ SINE=SIN(EPS) 00009080
+ DO 10 I=1,10 00009090
+ SINB=SIN(BETA(I)) 00009100
+ COSB=COS(BETA(I)) 00009110
+ SINL=SIN(LAMBDA(I)) 00009120
+ COSL=COS(LAMBDA(I)) 00009130
+ DELTA(I)=ASIN(SINB*COSE+COSB*SINE*SINL) 00009140
+ ALFA(I)=ATAN2(-SINB*SINE+COSB*COSE*SINL,COSB*COSL) 00009150
+C Lo conto antiorario da 0 a 360 00009160
+ IF(ALFA(I).LT.0) ALFA(I)=ALFA(I)+2*PIGREC 00009170
+ 10 CONTINUE 00009180
+ RETURN 00009190
+ END 00009200
+ SUBROUTINE ECLNOD(E,F,S,X) 00009210
+C ****************************************************************00009220
+C TRASFORMO DA COORDINATE ECLITTICHE RETTANGOLARI A COORDINATE 00009230
+C SUL PIANO ORBITALE, CON X LUNGO IL RAGGIO VETTORE 00009240
+C ****************************************************************00009250
+C QUESTA SUBROUTINE FA PARTE DEL PROGRAMMA DEGLI OROSCOPI 00009260
+C !!!!! LA PARTE PIU' NUMEROSA DEL COMITATO DI REDAZIONE ESPRIME 00009270
+C la sua perplessita' a riguardo della sottoindicata trasformazion00009280
+C di coordinate. Formulando un severo monito onde simili ghirigori00009290
+C non abbiano a ripetersi, declina ogni responsabilita' presente 00009300
+C e futura, non ritenendosi responsabile degli smarrimenti di 00009310
+C COORDINATE CHE SEGUIRANNO. !!!!!!!!! 00009320
+ IMPLICIT REAL*8 (A-H,O-Z) 00009330
+ DIMENSION A(3,3),X(3),S(3),E(6) 00009340
+ SINOF=SIN(E(2)+F) 00009350
+ COSOF=COS(E(2)+F) 00009360
+ SINOG=SIN(E(3)) 00009370
+ COSOG=COS(E(3)) 00009380
+ SINI=SIN(E(1)) 00009390
+ COSI=COS(E(1)) 00009400
+C TYPE*,'SINOF,COSOF,SINOG,COSOG,SINI,COSI' 00009410
+C TYPE*,SINOF,COSOF,SINOG,COSOG,SINI,COSI 00009420
+ A(1,1)=COSOF*COSOG-SINOF*SINOG*COSI 00009430
+ A(1,2)=COSOF*SINOG+SINOF*COSOG*COSI 00009440
+ A(1,3)=SINOF*SINI 00009450
+ A(2,1)=-SINOF*COSOG-COSOF*SINOG*COSI 00009460
+ A(2,2)=-SINOF*SINOG+COSOF*COSOG*COSI 00009470
+ A(2,3)=COSOF*SINI 00009480
+ A(3,1)=+SINOG*SINI 00009490
+ A(3,2)=-COSOG*SINI 00009500
+ A(3,3)=COSI 00009510
+C TYPE*,'A=',A 00009520
+C --------------------------------------------------------- 00009530
+ DO 10 I=1,3 00009540
+ S(I)=0. 00009550
+ DO 10 J=1,3 00009560
+ 10 S(I)=S(I)+A(I,J)*X(J) 00009570
+ RETURN 00009580
+ END 00009590
+ FUNCTION GIOR(ORA) 00009600
+C ******************************************************* 00009610
+C CONVERTE ORE IN FRAZIONI DI GIORNO 00009620
+C ! ORE ESPRESSE IN ORE,MINUTI 00009630
+C ******************************************************** 00009640
+ REAL*8 GIOR,ORA,A 00009650
+ A=INT(ORA) 00009660
+ GIOR=(A+(ORA-A)/60*100)/24 00009670
+ RETURN 00009680
+ END 00009690
+ FUNCTION GRARAD(G) 00009700
+C ******************************************* 00009780
+C SERVE PER CONVERTIRE I GRADI IN RADIANTI 00009710
+C ATTENZIONE! I GRADI DEVONO ESSERE COI DECIMALI 00009720
+C ******************************************* 00009780
+ REAL*8 GRARAD,G 00009730
+ GRARAD=G/180.D0*3.14159265358979323D0 00009740
+ RETURN 00009750
+ END 00009760
+ SUBROUTINE TRASF(X,A) 00009770
+C ******************************************* 00009780
+C TRASFORMA X COLLA MATRICE A 00009790
+C ************************************************ 00009800
+ IMPLICIT REAL*8 (A-H,O-Z) 00009810
+ DIMENSION Y(3),X(3),A(3,3) 00009820
+ DO 10 I=1,3 00009830
+ Y(I)=0.D0 00009840
+ DO 10 J=1,3 00009850
+ 10 Y(I)=A(I,J)*X(J)+Y(I) 00009860
+ GO TO 100 00009870
+ ENTRY INV(X,A) 00009880
+ DO 20 I=1,3 00009890
+ Y(I)=0.D0 00009900
+ DO 20 J=1,3 00009910
+ 20 Y(I)=A(J,I)*X(J)+Y(I) 00009920
+ 100 CONTINUE 00009930
+ DO 30 I=1,3 00009940
+ 30 X(I)=Y(I) 00009950
+ RETURN 00009960
+ END 00009970
+ SUBROUTINE JOD(DJ,G,M,JA) 00009980
+C ******************************************* 00009780
+C QUESTO E JODY,INVERSO DEL GIORNO GIULIANO COME SUBROUTINE 00010000
+C ******************************************* 00009780
+ REAL* 8 DJ,ALFA,A,B,G 00009990
+ IF(DJ.LE.0.) RETURN 00010010
+ DJ=DJ+0.5 00010020
+ IZ=INT(DJ) 00010030
+ IF(IZ-2299161) 3,5,5 00010040
+ 3 A=IZ 00010050
+ GO TO 6 00010060
+ 5 ALFA= INT((IZ-1867216.25)/36524.25) 00010070
+ A=IZ+1+ALFA-INT(ALFA/4.) 00010080
+ 6 CONTINUE 00010090
+ B=A+1524 00010100
+ IC=INT((B-122.1)/365.25) 00010110
+ ID=INT(365.25*IC) 00010120
+ IE=INT((B-ID)/30.6001) 00010130
+ F=DJ-IZ 00010140
+ G=B-ID-INT(30.6001*IE)+F 00010150
+ IF(IE.LT.13.5)M=IE-1 00010160
+ IF(IE.GT.13.5)M=IE-13 00010170
+ IF(M.LT.2.5) JA=IC-4715 00010180
+ IF(M.GT.2.5) JA=IC-4716 00010190
+ DJ=DJ-0.5 00010200
+C TYPE*,'JULIAN DAY= ',DJ 00010210
+C TYPE*,'GIORNO= ',G,' MESE= ',M,' ANNO= ',JA 00010220
+ RETURN 00010230
+ END 00010240
+ SUBROUTINE NODECL(E,XX,X,Y,Z) 00010250
+C ********************************************************** 00010260
+C MUTA VETTORE NELLE COORDINATE DEL PIANO DELL'ORBITA IN 00010270
+C COORDINATE SUL PIANO DELL'ECLITTICA 00010280
+C ********************************************************** 00010290
+C SUBROUTINE FACENTE PARTE DEL PROGRAMMA DEGLI OROSCOPI 00010300
+ IMPLICIT REAL*8(A-H,O-Z) 00010310
+ DIMENSION E(6),XX(3) 00010320
+C ---------------------------------------------------- 00010330
+ COSOG=COS(E(3)) 00010340
+ SINOG=SIN(E(3)) 00010350
+ COSI=COS(E(1)) 00010360
+ SINI=SIN(E(1)) 00010370
+C TYPE*,'NODECL:COSOG,SINOG,COSI,SINI' 00010380
+C TYPE*,COSOG,SINOG,COSI,SINI 00010390
+ XX(1)=X*COSOG-Y*SINOG*COSI+Z*SINOG*SINI 00010400
+ XX(2)=X*SINOG+Y*COSOG*COSI-COSOG*SINI*Z 00010410
+ XX(3)=Y*SINI+Z*COSI 00010420
+ RETURN 00010430
+ END 00010440
+ FUNCTION RADGRAD(R) 00010450
+C ******************************************* 00009780
+C SERVE PER CONVERTIRE I RADFIANTI IN GRADI 00010460
+C ATTENZIONE! FORNISCE I GRADI CON I DECIMALI 00010470
+C ******************************************* 00009780
+ REAL*8 RADGRAD,R 00010480
+ RADGRAD=R*180.D0/3.14159265358979323D0 00010490
+ RETURN 00010500
+ END 00010510
+ SUBROUTINE RETTECL(X,BETA,ALAMBD) 00010520
+C ****************************************** 00010530
+C DA COORDINATE RETTANGOLARI ECLITTICHE (X LUNGO GAMMA) 00010540
+C A BETA E LAMBDA , COORDINATE ECLITTICHE 00010550
+C **************************************************** 00010560
+ IMPLICIT REAL*8(A-H,O-Z) 00010570
+ DATA PIGREC/3.14159265358979323846/ 00010580
+ DIMENSION X(3,10),BETA(10),ALAMBD(10) 00010590
+ DO 10 I=1,10 00010600
+ ALAMBD(I)=ATAN2(X(2,I),X(1,I)) 00010610
+ BETA(I)=ASIN(X(3,I)/SQRT(X(1,I)**2+X(2,I)**2+X(3,I)**2)) 00010620
+C Lo conto antiorario dal gamma da 0 a 360 00010630
+ IF(ALAMBD(I).LT.0) ALAMBD(I)=ALAMBD(I)+2*PIGREC 00010640
+ 10 CONTINUE 00010650
+ RETURN 00010660
+ END 00010670
+ FUNCTION TESTG(Y,M,D) 00010680
+C ******************************************* 00009780
+C SERVE PER ELIMINARE LE DATE CHE NON ESISTONO 00010710
+C ******************************************* 00009780
+ IMPLICIT REAL*8 (A-H,O-Z) 00010690
+ INTEGER Y 00010700
+ DIMENSION MESE(12) 00010720
+ DATA MESE/31,29,31,30,31,30,31,31,30,31,30,31/ 00010730
+ IF(D.LT.0.OR.M.LT.0.OR.M.GT.12) GOTO 200 00010740
+ C=Y+M*1.E-2+D*1.E-4 00010750
+C GIORNI UCCISI DA GREGORIO 00010760
+ IF((C.GE.1582.1005).AND.(C.LE.1582.1014)) GOTO 200 00010770
+C MESE DI TROPPI DI' 00010780
+ IF(MESE(M).LT.D) GOTO 200 00010790
+C BISESTILE INESISTENTE 00010800
+ IF(INT(Y/4.)-Y/4..NE.0.AND.M.EQ.2.AND.D.GT.28) GOTO 200 00010810
+ IF(Y/400.-INT(Y/400.).NE.0.AND.M.EQ.2.AND.D.GT.28.AND. 00010820
+ 1 Y.GT.1582) GOTO 200 00010830
+ TESTG=0. 00010840
+ RETURN 00010850
+ 200 TESTG=1. 00010860
+ RETURN 00010870
+ END 00010880
+C 00010890
+
projects / frankie.git / commitdiff