Guida al FORTRAN applicato -...

199
Guida al FORTRAN applicato Lorenzo Zaninetti [email protected] Dipartimento di Fisica Generale Edizione III 10 febbraio 2018

Transcript of Guida al FORTRAN applicato -...

Guida al FORTRANapplicato

Lorenzo [email protected]

Dipartimento di Fisica GeneraleEdizione III

10 febbraio 2018

ii

Introduzione

Questa guida cerca di combinare due esigenze molto sentite al secon-do anno del corso di Laurea in Fisica: un corso di introduzione ad unlinguaggio di programmazione scientifica e alcuni esempi di programmi( funzionanti ! ) per l’ elaborazione dei dati raccolti durante le espe-rienze di laboratorio. Abbiamo quindi dedicato il capitolo uno ad unaintroduzione concisa al linguaggio di programmazione FORTRAN ed ilcapitolo due alla descrizione dei programmi per l’ elaborazione dei datiche possono essere usati dagli studenti in ambiente MS/DOS oppureLINUX oppure UNIX. Per soddisfare le esigenze degli studenti deglianni superiori riportiamo nel capitolo tre alcune subroutine di interessefisico-matematico che sono utilissime nei calcoli numerici di probabilitae statistica scritte in FORTRAN-77 standard. In realta questa guidapuo essere anche utile agli studenti degli anni superiori e quindi dedi-chiamo il capitolo quattro ai programmi del frattali e il capitolo cinquea quelli dei sistemi Monte Carlo. L’ appendice A riporta le funzioniintrinseche che possono essere richiamate dai programmi, l’ appendi-ce B alcune subroutines di interesse comune, l’ appendice C introducealcuni cenni del sistema operativo UNIX, l’ appendice D riporta i carat-teri ASCII, l’ appendice E i comandi del compilatore FORTRAN dellaMICROSOFT .

iii

CAPITOLO 0. INTRODUZIONE

iv

Indice

Introduzione iii

1 Il Fortran-77 11.1 Alcuni libri . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Caratteri e linee . . . . . . . . . . . . . . . . . . . . . . 21.3 Lunghezza della linea . . . . . . . . . . . . . . . . . . . 31.4 Linea di commento . . . . . . . . . . . . . . . . . . . . . 31.5 Linea iniziale e di continuazione . . . . . . . . . . . . . 31.6 Nomi simbolici . . . . . . . . . . . . . . . . . . . . . . . 51.7 Nomi locali e globali . . . . . . . . . . . . . . . . . . . . 51.8 Parole chiave . . . . . . . . . . . . . . . . . . . . . . . . 51.9 Etichette . . . . . . . . . . . . . . . . . . . . . . . . . . 61.10 Tipi di dati e costanti . . . . . . . . . . . . . . . . . . . 61.11 Unita di memoria . . . . . . . . . . . . . . . . . . . . . 61.12 Tipo INTEGER . . . . . . . . . . . . . . . . . . . . . . 61.13 Tipo REAL . . . . . . . . . . . . . . . . . . . . . . . . . 71.14 Tipo DOUBLE PRECISION . . . . . . . . . . . . . . . 71.15 Tipo COMPLEX . . . . . . . . . . . . . . . . . . . . . . 81.16 Tipo LOGICAL . . . . . . . . . . . . . . . . . . . . . . . 81.17 Tipo Character . . . . . . . . . . . . . . . . . . . . . . . 81.18 Vettori . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.19 Sottostringa CHARACTER . . . . . . . . . . . . . . . . 91.20 Espressioni . . . . . . . . . . . . . . . . . . . . . . . . . 91.21 Espressioni aritmetiche . . . . . . . . . . . . . . . . . . 101.22 Espressioni costanti intere . . . . . . . . . . . . . . . . . 111.23 Espressioni CHARACTER . . . . . . . . . . . . . . . . 111.24 Espressioni di relazione . . . . . . . . . . . . . . . . . . 121.25 Espressioni logiche . . . . . . . . . . . . . . . . . . . . . 131.26 Struttura del programma . . . . . . . . . . . . . . . . . 14

v

INDICE INDICE

1.27 Classi di istruzione . . . . . . . . . . . . . . . . . . . . . 151.28 Ordine delle istruzioni . . . . . . . . . . . . . . . . . . . 151.29 Sequenza di esecuzione . . . . . . . . . . . . . . . . . . . 171.30 Unita principale del programma . . . . . . . . . . . . . 171.31 Unita di programma FUNCTION . . . . . . . . . . . . 171.32 Funzioni di referenza . . . . . . . . . . . . . . . . . . . . 181.33 Funzioni istruzione . . . . . . . . . . . . . . . . . . . . . 191.34 Funzioni intrinseche . . . . . . . . . . . . . . . . . . . . 191.35 Sottoprogrammi SUBROUTINE . . . . . . . . . . . . . 191.36 Istruzione CALL . . . . . . . . . . . . . . . . . . . . . . 201.37 Argomenti della subroutine . . . . . . . . . . . . . . . . 201.38 Dimensione tramite PARAMETER . . . . . . . . . . . . 211.39 Dimensione tramite numeri . . . . . . . . . . . . . . . . 211.40 Dimensione tramite * . . . . . . . . . . . . . . . . . . . 221.41 Istruzione RETURN . . . . . . . . . . . . . . . . . . . . 221.42 Istruzione ENTRY . . . . . . . . . . . . . . . . . . . . . 221.43 Istruzioni di specificazione . . . . . . . . . . . . . . . . . 231.44 Istruzione IMPLICIT . . . . . . . . . . . . . . . . . . . 231.45 Istruzione PARAMETER . . . . . . . . . . . . . . . . . 231.46 Dichiarazioni dei vettori . . . . . . . . . . . . . . . . . . 241.47 Istruzioni tipologiche . . . . . . . . . . . . . . . . . . . . 251.48 Istruzione DIMENSION . . . . . . . . . . . . . . . . . . 251.49 Istruzione COMMON . . . . . . . . . . . . . . . . . . . 251.50 Istruzione EQUIVALENCE . . . . . . . . . . . . . . . . 261.51 Istruzione EXTERNAL . . . . . . . . . . . . . . . . . . 261.52 Istruzione INTRINSIC . . . . . . . . . . . . . . . . . . . 271.53 Istruzione SAVE . . . . . . . . . . . . . . . . . . . . . . 271.54 Istruzione DATA . . . . . . . . . . . . . . . . . . . . . . 281.55 Istruzioni di assegnazione . . . . . . . . . . . . . . . . . 281.56 Assegnazioni aritmetiche . . . . . . . . . . . . . . . . . . 281.57 Assegnazioni logiche . . . . . . . . . . . . . . . . . . . . 291.58 Assegnazioni CHARACTER . . . . . . . . . . . . . . . . 291.59 Istruzione ASSIGN . . . . . . . . . . . . . . . . . . . . . 291.60 Istruzioni di controllo . . . . . . . . . . . . . . . . . . . 301.61 Istruzione GO TO incondizionata . . . . . . . . . . . . . 301.62 Istruzione GO TO calcolata . . . . . . . . . . . . . . . . 301.63 Istruzione GO TO assegnata . . . . . . . . . . . . . . . 311.64 Blocco IF . . . . . . . . . . . . . . . . . . . . . . . . . . 311.65 Istruzione IF logica . . . . . . . . . . . . . . . . . . . . 32

vi

INDICE INDICE

1.66 Istruzione IF aritmetica . . . . . . . . . . . . . . . . . . 321.67 Istruzione DO . . . . . . . . . . . . . . . . . . . . . . . 321.68 Istruzione DO-WHILE . . . . . . . . . . . . . . . . . . . 331.69 Istruzione CONTINUE . . . . . . . . . . . . . . . . . . . 341.70 Istruzione STOP . . . . . . . . . . . . . . . . . . . . . . 351.71 Istruzione PAUSE . . . . . . . . . . . . . . . . . . . . . 351.72 Istruzioni di INPUT/OUTPUT . . . . . . . . . . . . . . 351.73 Nozioni generali . . . . . . . . . . . . . . . . . . . . . . 361.74 Record . . . . . . . . . . . . . . . . . . . . . . . . . . . 361.75 File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 371.76 Accesso ai file . . . . . . . . . . . . . . . . . . . . . . . . 371.77 File interni . . . . . . . . . . . . . . . . . . . . . . . . . 381.78 Istruzione OPEN . . . . . . . . . . . . . . . . . . . . . . 381.79 Istruzione CLOSE . . . . . . . . . . . . . . . . . . . . . . 391.80 Istruzione INQUIRE . . . . . . . . . . . . . . . . . . . . 401.81 Istruzioni READ/WRITE . . . . . . . . . . . . . . . . . 401.82 Istruzione REWIND . . . . . . . . . . . . . . . . . . . . 411.83 Istruzione BACKSPACE . . . . . . . . . . . . . . . . . . 421.84 Istruzione ENDFILE . . . . . . . . . . . . . . . . . . . . 421.85 La specificazione FORMAT . . . . . . . . . . . . . . . . 431.86 Descrittori ripetibili . . . . . . . . . . . . . . . . . . . . 431.87 Descrittori non ripetibili . . . . . . . . . . . . . . . . . . 45

2 Il Laboratorio 472.1 RETTA . . . . . . . . . . . . . . . . . . . . . . . . . . . 472.2 PARABO . . . . . . . . . . . . . . . . . . . . . . . . . . 502.3 PENDO . . . . . . . . . . . . . . . . . . . . . . . . . . . 522.4 CAVEN . . . . . . . . . . . . . . . . . . . . . . . . . . . 602.5 GAUSS . . . . . . . . . . . . . . . . . . . . . . . . . . . 742.6 TESTT . . . . . . . . . . . . . . . . . . . . . . . . . . . 792.7 TESTC2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

3 Probabilita e statistica 853.1 Gamma . . . . . . . . . . . . . . . . . . . . . . . . . . . 853.2 Fattoriale . . . . . . . . . . . . . . . . . . . . . . . . . . 863.3 Funzione beta . . . . . . . . . . . . . . . . . . . . . . . . 883.4 Gamma incompleta . . . . . . . . . . . . . . . . . . . . . 883.5 Funzione errori . . . . . . . . . . . . . . . . . . . . . . . 913.6 Distribuzioni . . . . . . . . . . . . . . . . . . . . . . . . . 92

vii

INDICE INDICE

3.6.1 Distribuzione del chiquadro . . . . . . . . . . . . 923.6.2 Distribuzione Beta incompleta . . . . . . . . . . . 933.6.3 Distribuzione di Student . . . . . . . . . . . . . . 95

3.7 Numeri random . . . . . . . . . . . . . . . . . . . . . . . 953.7.1 Generatori di sistema . . . . . . . . . . . . . . . . 953.7.2 Generatore Random - portatile . . . . . . . . . . 963.7.3 Generatore di interi . . . . . . . . . . . . . . . . . 973.7.4 Distribuzione esponenziale . . . . . . . . . . . . . 983.7.5 Distribuzione Gaussiana . . . . . . . . . . . . . . 99

3.8 Momenti . . . . . . . . . . . . . . . . . . . . . . . . . . . 1013.9 Test di Student . . . . . . . . . . . . . . . . . . . . . . . 1033.10 Fit con retta . . . . . . . . . . . . . . . . . . . . . . . . . 105

4 I frattali 1094.1 Henon mapping . . . . . . . . . . . . . . . . . . . . . . . 1094.2 Feigenbaum mapping . . . . . . . . . . . . . . . . . . . . 1124.3 Mandelbrot-set . . . . . . . . . . . . . . . . . . . . . . . 1154.4 Implementazione dei set . . . . . . . . . . . . . . . . . . 1164.5 Set di Julia e derivati . . . . . . . . . . . . . . . . . . . . 1204.6 Le trasformazioni affini . . . . . . . . . . . . . . . . . . . 1204.7 Il sistema L . . . . . . . . . . . . . . . . . . . . . . . . . 128

5 Sistemi Monte Carlo 1355.1 DLA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1355.2 Aggregazione ballistica . . . . . . . . . . . . . . . . . . . 1445.3 Percolazione . . . . . . . . . . . . . . . . . . . . . . . . . 1505.4 Catene di Voronoi . . . . . . . . . . . . . . . . . . . . . . 158

A Funzioni intrinseche 167

B Subroutine di interesse comune 175B.1 CLEANTOP . . . . . . . . . . . . . . . . . . . . . . . . 175B.2 MAXMIN . . . . . . . . . . . . . . . . . . . . . . . . . . 176

C Sistema operativo UNIX 177C.1 Alcuni comandi UNIX . . . . . . . . . . . . . . . . . . . 177C.2 Editor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178C.3 Compilazione+link/esecuzione . . . . . . . . . . . . . . . 179C.4 Shell Script . . . . . . . . . . . . . . . . . . . . . . . . . 179

viii

INDICE INDICE

C.5 Le librerie . . . . . . . . . . . . . . . . . . . . . . . . . . 180

D I caratteri ASCII 181

ix

INDICE INDICE

x

Capitolo 1

Il Fortran-77

Il linguaggio di programmazione FORTRAN, contrazione delle paroleFORmulaTRANslation, e un linguaggio per la soluzione di problemi computazio-nali. A causa della sua somiglianza con il familiare linguaggio aritmeticosemplifica enormemente la preparazione dei problemi da sottoporre allamacchina. Fu implementato per la prima volta nel 1956 e disponibileper gli utenti dal 1957. Negli anni tra il 1957 e il 1966 il linguaggio subıalcuni cambiamenti e fu standardizzato dall’ American National Stan-dards Institute (ANSI) nel 1966. Questa versione, essendo la quartadalla nascita del linguaggio, assunse il nome di FORTRAN IV. Il FOR-TRAN77 e l’ ultima implementazione ufficiale del FORTRAN-IV: com-missionato all’ ANSI per migliorare il linguaggio precedente e renderlopiu portabile, e stato finito nel 1977 ed e stato pubblicato come ANSIstandard X3-9-1978. Successivamente e stato adottato dalle maggioricompagnie e adesso funziona su tutti i PC e main frame. Infatti, inFORTRAN, sia i dati che le istruzioni sono organizzati in una sequen-za di istruzioni, chiamata normalmente ”codice sorgente”. Un codicesorgente scritto in FORTRAN puo essere eseguito su qualunque com-puter possegga un compilatore, con poche o nessuna modifica. Questeultime dipendono dal fatto che diverse case costruttrici forniscono im-plementazioni al linguaggio legate alle possibilita del tipo di computerche esse fabbricano. Occorre quindi, ad esempio, fare attenzione neltrasferimento di codici sorgente (o file) da UNIX/LINUX a personalcomputer.

1

1.1. ALCUNI LIBRI CAPITOLO 1. IL FORTRAN-77

1.1 Alcuni libri

Prima di continuare introduciamo brevemente una lista di libri chepossono essere utili per chi volesse avventurarsi in una trattazione piucompleta

1. Un’ introduzione ingegneristica alla programmazione stile ”Poli-tecnico” si trova in [Or 86]. Troverete anche una utile serie diesercizi e l’ associazione fra metalinguaggio e codice FORTRAN.

2. Una trattazione completa viene invece sviluppata in [El 85]. Ven-gono svolti o proposti molti esercizi.

3. Una trattazione sommaria viene invece svolta in [Do 83].

4. Un riassunto di tutte le istruzioni si trova invece in [Pa 83] : unmanualetto tascabile in inglese.

1.2 Caratteri e linee

L’ insieme dei caratteri in FORTRAN consiste di 26 lettere dalla A allaZ, le 10 cifre da 0 a 9 e dai seguenti 13 caratteri speciali:

+ ⇒ piu‘

− ⇒ meno

∗ ⇒ asterisco

⇒ spazio

= ⇒ eguale

/ ⇒ barra (slash)

( ⇒ parentesi sinistra

) ⇒ parentesi destra

. ⇒ punto decimale

, ⇒ virgola

: ⇒ due punti

′ ⇒ apostrofo

2

CAPITOLO 1. IL FORTRAN-77 1.3. LUNGHEZZA DELLA LINEA

< ⇒ minore

> ⇒ maggiore

$ ⇒ simbolo corrente

Attualmente il simbolo corrente puo avere un valore diverso a secondadella tastiera usata: non ha pero una funzione precisa in FORTRAN equindi si puo usare come simbolo di continuazione. Il carattere spazio eignorato dappertutto eccetto che nelle variabili CHARACTER e quindipuo essere usato a piacere per migliorare la leggibilita dei programmi.

1.3 Lunghezza della linea

Il compilatore fa uso solamente delle prime 72 colonne in ogni linea:dalla colonna 73 in poi si possono mettere dei simboli di identificazio-ne. Nelle linee di istruzione, le colonne da 1 a 5 sono riservate comeetichetta, la numero 6 ai segni di continuazione e quelle dal 7 al 72 perl’ istruzione stessa. Una linea di commento puo essere adoperata dallacolonna 2 in poi. Questa tipologia deriva dall’ impostazione originariatramite scheda, vedi figura 1.1.

1.4 Linea di commento

Qualsiasi linea con i caratteri ”C” o ”*” in prima colonna viene con-siderata una linea di commento. Le linee di commento non contanonulla ai fini della compilazione ed in esse si puo inserire qualsiasi ti-po di carattere: di solito vengono usate per annotazioni nelle varieparti del programma. La linea composta di spazi nelle sue prime 72colonne e equivalente ad uno spazio solo; ossia non conta ai fini dellacompilazione.

1.5 Linea iniziale e di continuazione

Una linea di istruzione e caratterizzata da uno spazio oppure da unozero nella colonna 6; dalla 1 alla 5 devono esserci degli spazi oppureuna etichetta. Qualsiasi istruzione, eccetto l’ END, puo continuaresu linee seguenti fino ad un totale di 20 linee ( una iniziale e 19 dicontinuazione ): le linee di commento possono essere frammezzate e

3

1.5. LINEA INIZIALE E DI CONTINUAZIONECAPITOLO 1. IL FORTRAN-77

Figura 1.1: La scheda

4

CAPITOLO 1. IL FORTRAN-77 1.6. NOMI SIMBOLICI

quindi non contano. Ogni linea di continuazione deve avere spazi dallacolonna 1 alla colonna 5 e il simbolo di continuazione in colonna 6.Il simbolo di continuazione puo essere qualsiasi carattere FORTRANeccetto, chiaramente, lo spazio o lo zero.

1.6 Nomi simbolici

Si puo usare un nome simbolico per identificare entita FORTRAN co-me costanti, variabili, vettori, funzioni e subroutines, blocchi commone unita di programma. Un nome simbolico e formato da una a 6 let-tere o numeri,con la condizione che il primo carattere sia una lettera.In alcuni casi, come per esempio in UNIX, i nomi possono essere piulunghi ma,di conseguenza, non sono trasportabili su tutti i sistemi. Laprima lettera determina il tipo di variabile ma questa regola puo essereignorata usando una particolare specificazione. Normalmente abbiamotipi INTEGER dalla lettera I alla N e REAL negli altri casi e questasituazione viene solitamente descritta dalla istruzione IMPLICIT.

1.7 Nomi locali e globali

I nomi di tutte le unita di programma, punti di entrata e blocchi COM-MON sono nomi globali, cioe comuni a tutto il programma: per questomotivo devono essere unici. Il nome di costanti, variabili, vettori e fun-zioni di istruzione sono locali, limitati all’ unita in uso, e quindi possonoessere adoperati per un’ altro scopo in un’ altra parte del programma.Sovente, pero, nome locale e globale devono essere differenti.

1.8 Parole chiave

Una parola chiave e una sequenza di caratteri con un significato pre-definito per il compilatore (per esempio istruzioni FORTRAN). Moltedi esse ( ad esempio DATA, IF, CLOSE, ERR ) hanno la forma di unnome simbolico valido: l’ interpretazione come parola chiave o comenome simbolico dipende solo dal contesto. Nel dubbio pero e meglioevitare i nomi delle istruzioni FORTRAN come nomi di variabili perevitare errori di compilazione o di esecuzione.

5

1.9. ETICHETTE CAPITOLO 1. IL FORTRAN-77

1.9 Etichette

Un’ istruzione di etichetta e un numero intero da 0 a 99999, senza segno,scritto nelle colonne da 1 a 5 della linea di istruzione. L’ etichetta nondeve essere allineata a destra e caratteri ” ” (spazio) e ”0” posti davantisono ignorati; deve essere unica nel programma e non puo essere messasu una linea di continuazione. Le etichette possono comparire non in or-dine numerico. L’ istruzione FORMAT deve sempre avere un’ etichetta;qualsiasi altra istruzione eseguibile puo avere un’ etichetta, compresal’ istruzione END.

1.10 Tipi di dati e costanti

Qualsiasi costante, variabile, vettore, espressione o funzione deve ap-partenere a uno di sei tipi di dati: INTEGER, REAL, DOUBLE PRE-CISION, COMPLEX, LOGICAL O CHARACTER. Il numero di tipipuo aumentare in alcune implementazioni: si ha, ad esempio, doppiaprecisione sulle variabili COMPLEX ( vedi UNIX-FORTRAN). Interi,reali, doppia precisione e complessi sono tipi aritmetici e la conversionetra di essi e fatta automaticamente in assegnazioni aritmetiche, DATA ePARAMETER. La conversione invece non e automatica nelle chiamateinterne tipo SUBROUTINE e FUNCTION.

1.11 Unita di memoria

Ogni variabile o elemento di vettore di tipo INTEGER, REAL, DOU-BLE PRECISION o LOGICAL occupa un’ unita di memoria numerica.Qualsiasi variabile o elemento di vettore in DOUBLE PRECISION oCOMPLEX occupa due unita di memoria numerica. Ogni variabileoppure elemento di vettore di tipo CHARACTER occupa un’ unita ditipo CHARACTER. Il rapporto fra le due unita dipende dalla macchinausata.

1.12 Tipo INTEGER

In una variabile INTEGER si memorizza un numero intero, che e rap-presentato esattamente nella macchina. Il range coperto dipende da

6

CAPITOLO 1. IL FORTRAN-77 1.13. TIPO REAL

macchina a macchina e nel UNIX-FORTRAN abbiamo:

integer ∗ 2 da − 32768 a 32768

integer ∗ 4 da − 2147483648 a 2147483648 .

In unita UNIX-FORTRAN una variabile integer*2 occupa 2 bytes men-tre quella integer*4 ne occupa 4.

1.13 Tipo REAL

Ogni variabile REAL memorizza un numero in virgola mobile, con osenza segno. Seppure lo zero e i numeri interi piccoli siano memorizzatiin maniera corretta, i numeri con virgola mobile sono approssimaticon una certa precisione che, seppure piccola, e finita e dipende dallamacchina. Nel UNIX-FORTRAN i numeri con virgola mobile variano:

real ∗ 4 da 0.79E − 38 a 1.78E + 38

real ∗ 8 da 0.56E − 308 a 0.9E + 308

Una costante reale deve contenere un punto ( che rappresenta lavirgola ) o un esponente che specifichi la potenza di 10 applicata allavariabile, o ambedue. Il segno e un optional da mettersi di fronte alnumero o all’ esponente. Esempi di costanti espresse in maniera valida:

−1. 3.14159 + 5E3 6.67E − 18

Sul UNIX-FORTRAN una variabile real *4 ( precisione singola ) occupa4 bytes.

1.14 Tipo DOUBLE PRECISION

La DOUBLE PRECISION e simile a quella REAL ma memorizza unnumero in virgola mobile con una precisione maggiore. Anche quılaprecisione dipende dalla macchina usata. Scriviamo qualche esempio:

−1D0 3.141592653589793238462643D0 + 5D3 6.67D − 18

Sul UNIX-FORTRAN una variabile real*8 (o doppia precisione) occupa8 bytes. Notiamo che l’ esponente viene preceduto dalla lettera D e devesempre essere indicato, anche se nullo.

7

1.15. TIPO COMPLEX CAPITOLO 1. IL FORTRAN-77

1.15 Tipo COMPLEX

Una variabile COMPLEX consiste in un paio ordinato di numeri reali,rappresentanti della parte reale e immaginaria di un numero complessoed e scritta con due costanti reali o intere separate da una virgolae contenute tra parentesi, la seconda delle quali rappresenta la parteimmaginaria:

(1, 50) (3.14,−2E3) (−.123,+4567)

1.16 Tipo LOGICAL

Una variabile LOGICAL memorizza solo due valori, .true. oppure .false.,ed occupa un byte: e quindi indicata per operazioni che richiedono moltibytes di memoria.

1.17 Tipo Character

I dati di tipo CHARACTER sono stati introdotti per la memorizza-zione e la manipolazione di testi. Ogni variabile CHARACTER puomemorizzare una stringa completa di caratteri, ma la sua occupazio-ne deve essere dichiarata all’ inizio del programma. In particolare inUNIX-FORTRAN bisogna specificare il numero di caratteri e quindi dibytes, per esempio

character ∗ 6 ⇒ ′ERRORE ′

character ∗ 15 ⇒ ′TROPPI ELEMENTI ′

Le stringhe sono definite da apostrofi (’) o virgolette (”). La costante diHollerith del FORTRAN-IV e stata completamente sostituita da questotipo di variabili ( per fortuna ! ).

1.18 Vettori

Un vettore e un set ordinato di dati (tutti dello stesso tipo) riferito adun singolo nome simbolico. Il numero di dimensioni di un vettore puoarrivare ad un massimo di sette. Nel FORTRAN77 qualsiasi dimensionepuo avere un limite inferiore diverso da uno e gli intervalli di variabilita

8

CAPITOLO 1. IL FORTRAN-77 1.19. SOTTOSTRINGA CHARACTER

devono essere dichiarati nelle istruzioni COMMON o DIMENSION.L’ elemento individuale viene chiamato specificando il nome simbolicoe un insieme di indici, uno per ogni dimensione.

Qualsiasi indice di un vettore deve avere un valore compreso fra idue limiti dichiarati in quella dimensione. L’ indice puo essere un interocomplesso quanto si vuole, puo includere la chiamata di un altro vettoreoppure chiamate a funzioni intrinseche od esterne. Riportiamo alcuniesempi:

matrice(0,−12)

hyper(i, j, l,m, n)

dlist(mod(npts, 4))

D(k(N(index)))

1.19 Sottostringa CHARACTER

Una sottostringa specifica un set di caratteri adiacenti dislocati in unavariabile CHARACTER. Una sottostringa consiste nella specificazio-ne dell’ elemento seguito da due interi separati da ”:” e rinchiusi fraparentesi. Facciamo un esempio concreto:

CHARACTER ∗ 8 regioni(1 : 23)

regioni(55) =′ PIEMONTE ′

regioni(55)(2 : 5) =′ IEMO′

1.20 Espressioni

Essendo un linguaggio essenzialmente computazionale, un gran nume-ro di istruzioni FORTRAN impiega espressioni. La valutazione di unaespressione da come risultato un singolo valore che puo essere usatoper definire una variabile, puo entrare in una decisione logica, puo es-sere scritto in un file, etc. La forma piu semplice di una espressionee un valore scalare: una costante o una variabile singola. Espressionipiu complesse vengono spesso formate con operazioni che coinvolgonouno o piu operandi. Ci sono quattro tipi di espressioni disponibili inFORTRAN 77: aritmetiche, CHARACTER, di relazione e logiche.

9

1.21. ESPRESSIONI ARITMETICHE CAPITOLO 1. IL FORTRAN-77

1.21 Espressioni aritmetiche

Un’ espressione aritmetica e usata per specificare un calcolo numericoche deve essere svolto dal calcolatore. Consiste di uno o piu operan-di ( interi, reali, doppia precisione oppure complessi ) combinati conoperatori aritmetici. Un’ espressione puo anche contenere una sottoe-spressione in parentesi e referenze a funzioni standard. Il risultato diun’ espressione aritmetica e un valore numerico del tipo specificato.Elenchiamo le possibili operazioni:

somma(identita) ⇒ +

sottrazione(negazione) ⇒ −

moltiplicazione ⇒ ∗

divisione ⇒ /

elevazione a potenza⇒ ∗∗

L’ ordine di valutazione di un’ espressione e:

1. sottoespressioni in parentesi,

2. funzioni standard,

3. esponente,

4. moltiplicazione o divisioni,

5. addizioni, sottrazioni o negazioni.

In ogni gruppo la valutazione procede da sinistra verso destra eccet-to che per gli esponenti che sono valutati da destra a sinistra. CosiccheA/B/C e equivalente ad (A/B)/C mentre A**B**C e equivalente adA** (B**C). Si possono mescolare tutte le combinazioni dei quattro ti-pi fondamentali tra di loro eccetto che quelle COMPLEX con quelle diDOUBLE PRECISION. Se gli operandi di una espressione aritmeticasono dello stesso tipo il risultato sara anch’ esso dello stesso tipo. Se itipi sono diversi, viene applicata ad uno di essi un tipo di conversioneimplicita. La direzione adoperata sara:

INTEGER⇒ REAL⇒ COMPLEX o DOUBLE PRECISION

10

CAPITOLO 1. IL FORTRAN-77 1.22. ESPRESSIONI COSTANTI INTERE

La sola eccezione a queste regole di conversione e l’ elevazione a potenzacon esponente intero. Se la potenza e negativa, si valuta l’ inverso dellapotenza positiva:

2 ∗ ∗(−3) ⇒ 1/(2.0 ∗ ∗3) ⇒ 1/8.0⇒ 0.125

Altrimenti la potenza e calcolata con i logaritmi cosicche A**B saracalcolato con la formula EXP(B * LOG (A)) dove EXP e LOG sonofunzioni standard del compilatore. La divisione fra interi produce untroncamento della parte frazionale. Esempi:

−8/3 ⇒ −2 , 2 ∗ ∗(−3) ⇒ 1/(2 ∗ ∗3) ⇒ 1/8 ⇒ 0

1.22 Espressioni costanti intere

Le espressioni aritmetiche in cui tutti gli operandi siano costanti o nomisimbolici, vengono chiamate espressioni aritmetiche costanti. Un’ espres-sione costante intera puo essere usata per dichiarare i limiti di un vettoreo la lunghezza di una variabile CHARACTER o usata per fissare i valoridi una costante attraverso l’ istruzione PARAMETER. Per esempio:

PARAMETER (NPT = 50, NTOTAL = NPT ∗ ∗2)

COMPLEX COEFFS (NPTS ∗ (NPTS − 1))

ARRAY (NTOTAL(2))

1.23 Espressioni CHARACTER

Una espressione CHARACTER produce un risultato di tipo CHARAC-TER ed e ottenuta usando operatori e operandi di tipo CHARACTER.Gli operandi di tipo CHARACTER possono essere: un valore scalare,un elemento di vettore, una sottostringa, una espressione CHARAC-TER chiusa tra parentesi, il risultato di una funzione CHARACTER.Esiste un unico operando CHARACTER denotato con // che permettedi concatenare due stringhe. Esempio:

A = ′ESEMPIO DI ′

B = ′PROGRAMMA′

11

1.24. ESPRESSIONI DI RELAZIONE CAPITOLO 1. IL FORTRAN-77

C = ′SUBROUTINE ′

A//B = ′ESEMPIO DI PROGRAMMA′

A//C = ′ESEMPIO DI SUBROUTINE ′

Ambedue gli apostrofi sono richiesti .Il valore della variabile CHARACTER e la stringa di caratteri fra

gli apostrofi. Il suddetto valore non include gli apostrofi, include invecetutti gli spazi ed eventuali tabs inseriti fra i due apostrofi . Nellavariable CHARACTER il carattere apostrofo e rappresentato da dueapostrofi consecutivi con nessun spazio o carattere fra di loro (”) . Lalunghezza della variabile CHARACTER e il numero di caratteri fragli apostrofi, eccetto che due apostrofi consecutivi rappresentano unsingolo apostrofo . La lunghezza di una variabile CHARACTER deveessere nel range da 1 a 2000 .Esempi

I seguenti esempi rappresentano delle variabili CHARACTER valideed invalide .Valide’Cosa hai detto ?’’L”abbiamo trovato’’La risposta media e stata ”12:00” ’

Invalide Spiegazione’La formula manca l’ apostrofo finale” la variabile CHARACTER deve

contenere almeno un carattere

1.24 Espressioni di relazione

Un’ espressione di relazione da un risultato di tipo logico (vero o falso)ed e formata da espressioni aritmetiche o CHARACTER e operatori direlazione. Abbiamo sei tipi di operatori di relazione:

.LT. ⇒ minore di

.LE. ⇒ minore di o uguale a

12

CAPITOLO 1. IL FORTRAN-77 1.25. ESPRESSIONI LOGICHE

.EQ. ⇒ uguale a

..GT. ⇒ maggiore di

.GE. ⇒ maggiore di o uguale a

.NE. ⇒ non uguale a

Le variabili confrontate devono essere omogenee e ovviamente in campocomplesso potremo solamente usare l’espressione .eq. o.ne..

Inoltre nei test di uguaglianza su variabili REAL o DOUBLE PRE-CISION, bisogna tenere conto della precisione adoperata. Il confrontotra variabili CHARACTER avviene carattere per carattere secondo lasequenza ASCII. Dato che la lettera ”A” precede la lettera ”B” nelcodice ASCII, ”A” e minore di ”B”. Inoltre le lettere maiuscole hannovalori minori delle minuscole.

1.25 Espressioni logiche

Un’ espressione logica consiste di operandi e operatori logici e ritornaun valore . true./.false .. Essa puo essere usata nell’ espressione IFo blocco IF. Dato il termine logico lterm possiamo avere le seguentiespressioni logiche:

(lterm)

(.not. lterm)

(lexp lop lterm)

dove LTERM e la costante logica, variabile oppure elemento di vet-tore, LEXP un’ espressione logica e LOP uno degli operatori logi-ci.AND.,.OR.,.EQV., o.NEQV.. Riportiamo l’ ordine di valutazionedi un’ espressione di tipo logico:

1. Espressione aritmetica

2. Espressioni di relazione

3. operatore.NOT.

4. operatore.AND.

5. operatore.OR.

13

1.26. STRUTTURA DEL PROGRAMMA CAPITOLO 1. IL FORTRAN-77

Tabella 1.1: Tabella sugli operatori OR e AND

L1 ed L2 sono espressioni logiche.OR. e.AND. sono operatori logici

L1 L2 L1.OR.L2 L1.AND.L2

.true. .true. .true. .true.

.true. .false. .true. .false..false. .true. .true. .false..false. .false. .false. .false.

6. operatori.EQV. e.NEQV.

L’ effetto degli operatori.OR. da come risultato .true. se e vero unodei due operandi mentre.AND. da come risultato.true. se sono verificatitutti e due, vedi per esempio la seguente tabella:

L’ operatore.EQV. da come risultato.true. se entrambi gli operan-di hanno lo stesso valore ( cioe entrambi veri oppure entrambi falsi )mentre il secondo.NEQV. da il risultato.true. se hanno valore opposto,come indicato nella tabella che segue:

Attenzione che l’ operatore.EQV. (che esprime equivalenza ) logicae.NEQV. (non - equivalenza) non possono essere adoperati nella stes-sa espressione. Le parentesi vanno sempre utilizzate per semplificarel’ ordine di espressioni complicate.

1.26 Struttura del programma

Un programma deve contenere un corpo principale, detto main. Puopoi contenere diverse sottounita che possono essere di tre tipi: SUB-ROUTINE, FUNCTION e blocco DATA (procedura non eseguibile).

14

CAPITOLO 1. IL FORTRAN-77 1.27. CLASSI DI ISTRUZIONE

Tabella 1.2: Tabella sugli operatori EQV e NEQV

L1 ed L2 sono espressioni logicheEQV. e.NEQV. sono operatori logici

L1 L2 L1.EQV.L2 L1.NEQV.L2

.true. .true. .true. .false.

.true. .false. .false. .true..false. .true. .false. .true..false. .false. .true. .false.

1.27 Classi di istruzione

Possiamo chiamare istruzioni tipologiche INTEGER, REAL, DOUBLEPRECISION,COMPLEX, LOGICAL e CHARACTER mentre IMPLICIT, PARA-METER, DIMENSION, COMMON, EQUIVALENCE, EXTERNAL,SAVE sono classificate come istruzioni di specificazione. PROGRAM,SUBROUTINE, FUNCTION, blocco DATA, FORMAT sono inveceistruzioni di funzione e non sono eseguibili. Le istruzioni di assegnazio-ne ( incluso ASSIGN ), tutte le forme di IF e GO TO, e inoltre ELSE,ELSEIF, END IF,DO, CONTINUE, CALL, RETURN, READ, WRI-TE, PRINT, OPEN, CLOSE, INQUIRE, REWIND, BACKSPACE,END FILE, STOP, PAUSE ed END sono classificate come istruzionieseguibili.

1.28 Ordine delle istruzioni

Nella scrittura del codice sorgente devono essere rispettate alcune regoledi priorita con un certo ordine, vedi figura 1.2 .♣ PROGRAM NOME♣ Istruzioni di specificazione♣ ........................♣ Istruzioni eseguibili♣ .......................

15

1.28. ORDINE DELLE ISTRUZIONI CAPITOLO 1. IL FORTRAN-77

+-------+--------------------------------------------------------+

| | Istruzione OPTIONS |

| |--------------------------------------------------------|

| | Istruzione PROGRAM, FUNCTION, SUBROUTINE,o BLOCCO DATA |

| |--------+-----------------------------------------------|

|COMMENT| | Istruzione IMPLICIT NONE |

| Linee,| |-------------------------------+--------------|

|INCLUDE|NAMELIST,| Istruzione IMPLICIT | |

| Istru-| FORMAT, |------+------------------------| |

| zioni,| & | | Altre Istruzioni | ISTRUZIONE |

|& Gen- | ENTRY | DATA | di specificazione | PARAMETER |

| erali | Istru- |Istru-| Istruzione DICTIONARY | |

|Diret- | zioni | zioni|------------------------+--------------|

| tive | | | Istruzioni Function Definizioni |

| | | |---------------------------------------|

| | | | Istruzioni Eseguibili |

|-------+---------+------+---------------------------------------|

| Istruzione END |

+----------------------------------------------------------------+

Figura 1.2: Le priorita nelle istruzioni

♣ STOP♣ END

Questo significa, per esempio, che se si usa una istruzione di tipoIMPLICIT,questa deve precedere tutte le altre eccetto che PARAME-TER, inseribile ovunque. Le linee di commento possono essere inseri-te in ogni punto eccetto che dopo l’ istruzione END, che deve esseresempre l’ ultima del programma. Nella figura che precede le linee ver-ticali separano le istruzioni che possono essere mescolate. Per esempiol’ istruzione DATA puo essere mescolata con istruzioni eseguibili. Lelinee orizzontali invece separano istruzioni che non possono essere me-scolate; per esempio le istruzioni di dichiarazione non possono esseremescolate con istruzioni eseguibili.

16

CAPITOLO 1. IL FORTRAN-77 1.29. SEQUENZA DI ESECUZIONE

1.29 Sequenza di esecuzione

L’ esecuzione parte con la prima istruzione eseguibile del programmamain e continua sequenzialmente, soggetta pero alle istruzioni di con-trollo come GO TO, IF, CALL, etc.. L’ esecuzione finisce con l’ istruzio-ne END del programma principale o con uno STOP in qualsiasi partedel programma. Certi errori tipo OVERFLOW possono creare una fineprematura ma segnalata mentre altri, tipo superamento della memoriadisponibile, provocano una fine non segnalata.

1.30 Unita principale del programma

Il main incomincia con l’ istruzione PROGRAM e puo contenere qual-siasi altraistruzione eccetto che: SUBROUTINE, FUNCTION oppure RETURN.La forma dell’ istruzione e:♣ PROGRAM pnome

dove pnome e un nome simbolico. Questa istruzione non e adoperatain altre parti del programma e quindi ha un significato principalmentedi documentazione.

1.31 Unita di programma FUNCTION

Una funzione ritorna un valore al punto, nell’ ambito di una espres-sione, in cui e stata chiamata. Quest’ unita comincia con l’ istruzio-ne FUNCTION e puo contenere qualsiasi altra istruzione eccetto chePROGRAM, SUBROUTINE o blocco DATA. Abbiamo due possibiltadi definizione:♣ FUNCTION fnome ( a1, a2,.... an )

oppure♣ tipo FUNCTION fnome (a1, a2,... an )

dove tipo va scelto fra INTEGER, REAL, DOUBLE PRECISION,COMPLEX, LOGICAL, CHARACTER o CHARACTER * len; fnomee il nome simbolico della funzione e a1, a2,.... an sono gli argomentitrasmessi. Si puo anche non trasferire alcun argomento ma le parentesidevono ancora essere presenti. E meglio evitare che nella FUNCTIONci siano operazioni di PRINT/WRITE perche non funzionerebbero ed

17

1.32. FUNZIONI DI REFERENZA CAPITOLO 1. IL FORTRAN-77

inoltre bisogna essere estremamente prudenti nel modificare gli argo-menti indicati nell’ istruzione COMMON. Come esempio riportiamo lafunzione MYSIN che calcola il seno del primo argomento con lo svi-luppo in serie di Taylor e ritorna il numero dei termini dello svilupponecessari ad avere una precisione di 10−7.♣ REAL FUNCTION MYSIN(P1,P2)♣ REAL P1♣ DOUBLE PRECISION FACTRL♣ INTEGER P2♣ MYSIN = P1♣ DO 10 P2 = 1,100♣ DELTA = (-1)**P2 * P1**(2*P2+1) / FACTRL(2*P2+1)♣ IF (ABS(DELTA).LT. 1E-7) RETURN♣10 MYSIN = MYSIN + DELTA♣ RETURN♣ END

Bisogna fare attenzione a dichiarare la funzione chiamata nel pro-gramma principale; nel nostro caso dovremo mettere♣ REAL MYSIN

se poi la funzione e inserita in qualche libreria dovremo mettere♣ EXTERNAL MYSIN

1.32 Funzioni di referenza

Una funzione di referenza puo occorrere come operando in qualsiasiespressione appropriata ed ha la forma del tipo:♣ fnome ( a1,....., an)

dove fnome e il nome della funzione. Quella piu usate sono le fun-zioni matematiche che variano fortemente da macchina a macchina. Ri-portiamo l’ esempio dei numeri random fra 0 ed 1 in UNIX-FORTRAN:♣ X = RAN (IGEN)

dove IGEN deve essere un intero grande meglio se dispari. Le funzio-ni intrinsiche non vanno dichiarate perche sono inserite nel compilatore.Per maggiori dettagli vedere l’ appendice A sulle funzioni intrinsiche.

18

CAPITOLO 1. IL FORTRAN-77 1.33. FUNZIONI ISTRUZIONE

1.33 Funzioni istruzione

Una istruzione di funzione puo essere usata per specificare una relazionesemplice ed e usata in una sola unita. Ovviamente essa deve precederequalsiasi istruzione di esecuzione. E definita tramite una via simile alleistruzioni di assegnazione:♣ fnome (a1,...,an) = exp

dove fnome e il nome simbolico della funzione ed e seguito dagliargomenti fittizi, exp e una espressione che viene valutata utilizzandogli argomenti fittizi. Questi ultimi sono locali e possono essere usatiindipendentemente nel resto del programma. Esempi:♣ LOGICAL DIGIT♣ CHARACTER CHAR♣ DIGIT (CHAR ) = CHAR.GE.’0’.AND.CHAR.LE.’9’♣ VAB (X,Y,Z) = SQRT (XX**2 +Y**2 +Z**2)

1.34 Funzioni intrinseche

Il FORTRAN contiene una quantita considerevole di funzioni incor-porate, note come funzioni intrinseche, che vengono automaticamentemesse a disposizione in qualsiasi parte del programma esse venganochiamate. La maggior parte di queste funzioni ha un nome genericoche elimina la necessita di cambiare nome in base al tipo di parametrousato. Su varie macchine troviamo poi funzioni intrinseche diverse aseconda del tipo di variabile utilizzata. Visto il loro numero notevolerimandiamo il lettore all’ appendice A per una lista aggiornata.

1.35 Sottoprogrammi SUBROUTINE

Una subroutine e una procedura separata, definita come esterna rispet-to all’ unita di programma che la chiama ed e specificata in una unitasottoprogramma. L’ istruzione SUBROUTINE seguita o meno dagliargomenti e l’ istruzione iniziale e in seguito si possono inserire tuttele istruzioni eccetto che PROGRAM, FUNCTION o blocco DATA. Laforma puo essere:♣ SUBROUTINE snome (a1,....,an)

oppure♣ SUBROUTINE snome

19

1.36. ISTRUZIONE CALL CAPITOLO 1. IL FORTRAN-77

dove snome e il nome simbolico della subroutine e a1,....., an sonogli argomenti di entrata / uscita. La SUBROUTINE viene chiamatacon l’ istruzione CALL e ovviamente al suo interno non puo chiamarese stessa direttamente o indirettamente. Le informazioni si possonotrasferire al di fuori della subroutine mediante i suoi argomenti, i blocchiCOMMON oppure file esterni. Si ritorna al main mediante l’ istruzioneEND oppure RETURN.

1.36 Istruzione CALL

Questa istruzione e usata per eseguire una procedura specificata da unasubroutine. Le forme sono del tipo:♣ CALL SUBROUTINE snome (a1,....,an)

oppure♣ CALL SUBROUTINE snome

1.37 Argomenti della subroutine

Normalmente l’ istruzione iniziale SUBROUTINE contiene una listadi parametri formali chiusi fra parentesi. La chiamata CALL ad unasubroutine contiene una lista di parametri attuali che devono concor-dare esattamente in numero e tipo con i parametri formali. Durantel’ esecuzione di un’ istruzione CALL, le locazioni di memoria dei para-metri attuali vengono passate alla subroutine in modo tale da permet-tere ai parametri formali di riferirsi alle stesse locazioni di memoria deiparametri attuali. Questi ultimi possono essere costituiti da:

1. una costante

2. un nome di variabile

3. il nome di un elemento di un vettore o del vettore stesso

4. una sottostringa di caratteri

5. un’ espressione

6. il nome una funzione intrinseca

7. il nome di una procedura esterna ( subroutine o funzione )

20

CAPITOLO 1. IL FORTRAN-771.38. DIMENSIONE TRAMITE PARAMETER

8. il nome di una procedura formale

9. uno specificatore di ritorno alternativo ( vedi sotto )

I parametri formali possono essere invece costituiti da:

• a) il nome di una variabile ( per i tipi 1 - 5 )

• b) il nome di un vettore ( per il tipo 3 )

• c) il nome di una procedura formale ( per i tipi 7 e 9 )

• d) un asterisco ( per il tipo 10 )

Se il parametro formale e una variabile di tipo CHARACTER, deveessere naturalmente dichiarato come tale e gli deve essere attribuitauna certa lunghezza che deve essere inferiore o uguale alla lunghezzadel parametro attuale; se e inferiore, verranno trattati come parametroformale solo i ’lun’ ( lun = lunghezza dichiarata nella SUBROUTINE) caratteri piu a sinistra del parametro attuale. Particolare attenzioneva posta alle dimensioni dei vettori e possiamo avere tre tipi diversi diconfigurazioni:

1.38 Dimensione tramite PARAMETER

♣ SUBROUTINE LAVORO (VEC,N)♣ DIMENSION VEC (1:N)

dove VEC e il vettore di dimensione da uno ad N e N e definito nelmain tramite l’ istruzione PARAMETER.

1.39 Dimensione tramite numeri

♣ SUBROUTINE LAVORO (VEC)♣ DIMENSION VEC (1:1000)

dove VEC e il vettore di dimensione da uno ad a 1000 e nel MAIN edefinito nella stessa maniera.

21

1.40. DIMENSIONE TRAMITE * CAPITOLO 1. IL FORTRAN-77

1.40 Dimensione tramite *

♣ SUBROUTINE LAVORO (A,B,C)♣ DIMENSION A (1:*),B(0:*),C(-4:*)Questo tipo di dichiarazione e detto ”dichiarazione di vettore a di-mensioni presunte”, in quanto si suppone che il vettore sia abbastanzagrande per tutti i riferimenti ad esso fatti nel corso del programma enon ha dimensioni definite. Se invece abbiamo l’ istruzione in bloc-chi COMMON abbiamo un associazione fra elementi corrispondenti; siabbia ad esempio♣ SUBROUTINE SUB ( X,Y,Z)♣ COMMON C

e nell’ unita principale si abbiano le istruzioni seguenti:♣ COMMON S♣ CALL SUB (X1,Y1,Z1 )

Cosicche sia X,Y,Z sono associati con X1, Y1, Z1 mentre C e asso-ciata con S.

1.41 Istruzione RETURN

Questa istruzione fa terminare l’ esecuzione di una SUBROUTINE odi una FUNCTION e restituisce il controllo all’ unita di chiamata. Laforma normale e:♣ RETURN

Puo anche essere chiamata all’ interno di un if /end if.

1.42 Istruzione ENTRY

Qualche volta puo essere utile combinare due o piu subroutines o fun-zioni in una sola per potere trarre vantaggio da una grossa porzione dicodice in comune. Essa assume la forma del tipo:♣ ENTRY enome (a1,....,an)

dove enome e il nome di entrata seguito dai relativi argomenti, ilcui numero deve essere uguale al numero di argomenti nella subroutineo funzione corrispondente. L’ uso di punti alternativi di ingresso e inconflitto con i principi di programmazione modulare ma puo risultareconveniente in alcuni casi particolari. Bisogna evitare comunque chel’ istruzione ENTRY compaia in una qualsiasi struttura a blocchi

22

CAPITOLO 1. IL FORTRAN-77 1.43. ISTRUZIONI DI SPECIFICAZIONE

1.43 Istruzioni di specificazione

Le istruzioni di specificazione sono usate per definire le proprieta dientita simboliche, costanti, variabili, vettori, etc., usate in una unita diprogramma. Per questo motivo vengono anche chiamate istruzioni didichiarazione e sono raggruppate nella sezione iniziale del codice sorgen-te, dopo l’ istruzione PROGRAM e prima di tutte le altre. Ricordiamoche non sono istruzioni eseguibili. In questo capitolo si esamineran-no anche le istruzioni DATA, utilizzate per stabilire i valori iniziali divariabili e vettori presenti in un programma FORTRAN.

1.44 Istruzione IMPLICIT

Il tipo di dati delle variabili o vettori usati nel programma puo esserestabilito tramite un’ istruzione tipologica. La regola standard e chei nomi simbolici con la lettera iniziale da I ad N specificano variabiliintere mentre gli altri sono reali. Con l’ istruzione IMPLICIT si puomodificare questa situazione. Questo puo essere utile per poter utiliz-zare allo stesso tempo variabili in precisione semplice e in doppia nellestesse unita. Se usiamo il tipo CHARACTER deve essere specificataanche la lunghezza con un *len. Esempi:♣ IMPLICIT REAL *4 (A-H,O-Z)

♣ oppure♣ IMPLICIT REAL *8 (A-H,O-Z)

♣ oppure♣ IMPLICIT REAL *4 (C-H,O-Z)♣ IMPLICIT REAL *8 (A)♣ IMPLICIT CHARACTER * 80 ( B)

In questo ultimo esempio tutte le variabili saranno in semplice pre-cisione eccetto che quelle inizianti con A ( precisione doppia ) e B (variabile CHARACTER di lunghezza 80 ).

1.45 Istruzione PARAMETER

Con questa istruzione si assegna un nome logico ad una costante, chepuo essere un numero oppure un’ espressione contenente altre costan-ti. Si possono trasportare le variabili cosı definite in altre unita senza

23

1.46. DICHIARAZIONI DEI VETTORI CAPITOLO 1. IL FORTRAN-77

modificarle: il compilatore segnalerebbe in questo caso qualcosa di ano-malo. Il tipo di costante o espressione deve chiaramente essere lo stessodel nome simbolico. Se il nome di un parametro CHARACTER hala lunghezza specificata con un *(*), vale la lunghezza dichiarata inPARAMETER. Esempio:♣ PARAMETER (LOGN=8, NPT = 2** LOGN)♣ DIMENSION INDEX (1:LOGN), X(0:NPT)♣ CHARACTER *(*) TITLE♣ PARAMETER ( TITLE =’Tabella dei risultati’)♣ PARAMETER ( PI = 3.14159,PI2 = 2*PI )

Esse risultano molto utili per dimensionare i vettori.

1.46 Dichiarazioni dei vettori

Il nome di un vettore puo comparire in piu di una istruzione di speci-ficazione ma la sua dimensione deve essere dichiarata solo una volta.Questo puo essere fatto in una istruzione tipologica, in un COMMONoppure in un DIMENSION. I vettori possono avere fino a sette dimen-sioni, per ognuna delle quali si deve specificare un limite inferiore esuperiore. Se non e specificato diversamente, il limite inferiore e uno.Esempi:♣ PARAMETER (NPUNTI = 1000 )♣ REAL 8 4 (NPUNTI),Y(-5:5,0:NPUNTI -1),DUE(2,2,2,2)♣ LOGICAL * 1 CUBO♣ DIMENSION CUBO (1:NPUNTI,1:NPUNTI,1:NPUNTI)

Se invece il vettore si trova in una subroutine abbiamo due possibi-lita: dimensione regolabile oppure fissa. Nel primo caso le dimensionisono regolate dalla lista degli argomenti della subroutine oppure tramiteil blocco COMMON:♣ SUBROUTINE CALCOLA (X,NX,Y,NY, R,F)♣ IMPLICIT REAL *8 (A-H,O-Z)♣ COMMON NR,NF♣ DIMENSION X (1:NX),Y(1:NY),F(0:NF),R(0:NR)

Se invece vogliamo usare un vettore di dimensione fissa possiamointrodurre un ”*” che rappresenta una dimensione indefinita. Esempio:♣ FUNCTION ESPANDE (XX,YY)♣ IMPLICIT REAL *8 XX, YY♣ DIMENSION XX (*),YY(0:*)

24

CAPITOLO 1. IL FORTRAN-77 1.47. ISTRUZIONI TIPOLOGICHE

1.47 Istruzioni tipologiche

Queste sono usate per alterare ( o confermare ) il tipo di dati associaticon un nome simbolico. Possono anche essere usate per specificare ledimensioni dei vettori. Una volta associato un tipo ad un nome sim-bolico, questo rimane invariato all’ interno dell’ unita di programmain cui e stato definito. Possiamo avere tipi LOGICAL, CHARACTERe aritmetici. Usando l’ istruzione CHARACTER bisogna ovviamen-te sempre specificare la lunghezza. Nel caso dell’ argomento fittiziodi una SUBROUTINE possiamo mettere un *(*) che rappresenta unalunghezza e dimensione variabili. Esempio:♣ DOUBLE PRECISION X, FUNZ, Z(100)♣ CHARACTER * 30 NOME, INDIRIZZI (5)*20♣ LOGICAL * 1 CUBO (1:100,1:100)

1.48 Istruzione DIMENSION

E usata per dare un nome simbolico e specificare la dimensione di unvettore. Il tipo di dati e implicito nel nome usato. La forma e del tipo:♣ DIMENSION A1(D1),A2(D2).....

dove ogni An e un vettore e Dn specifica la dimensione; a sua voltaDn puo essere scomposto in limite superiore e inferiore.

1.49 Istruzione COMMON

Le istruzioni COMMON sono usate per dichiarare blocchi di variabiliin comune con due o piu unita del programma, secondo un ordine dimemorizzazione consistente. Questo e possibile attraverso una struttu-ra di dati FORTRAN chiamata ” blocco COMMON”, che costituisceun blocco di memoria contiguo e assume la forma:♣ COMMON /NOME/ nlist

dove ”nome” e il nome globale del blocco e ”nlist” e la lista deinomi delle variabili, dei vettori e dei dichiaratori di vettori locali. Sinota che un blocco COMMON non possiede un tipo. Una volta stabilitol’ ordine di memorizzazione in un blocco COMMON, qualsiasi unita diprogramma che dichiari lo stesso blocco COMMON, puo riferirsi ai datiivi memorizzati senza dover passare nomi simbolici attraverso liste diargomenti. L’ istruzione COMMON e un istruzione di specificazione e

25

1.50. ISTRUZIONE EQUIVALENCE CAPITOLO 1. IL FORTRAN-77

deve venire quindi dopo tutte le istruzioni IMPLICIT e precedere tuttequelle DATA. Esempio:♣ INTEGER *2 ETA (1:50)♣ REAL *4 LIVELLO (1:50,1:6)♣ COMMON /ESAM/NUM,LIVELLO,ETA

In questo caso nel blocco COMMON /ESAM/ abbiamo 1 + 300 +50 = 351 unita.

1.50 Istruzione EQUIVALENCE

L’ istruzione EQUIVALENCE fornisce un mezzo affinche due o piuvariabili possano occupare la stessa dislocazione in memoria. Assumela forma:♣ EQUIVALENCE (nlista1),(nlista2)

Questo significa che la memorizzazione per tutti gli elementi del-la lista avviene con la stessa locazione di partenza. Gli elementi inquestione possono essere di tipo e lunghezza differente. Facciamo unesempio:♣ REAL *4 R(2), RL,IM♣ COMPLEX *8 C♣ EQUIVALENCE ( R,RL,C),(R(2),IM)

Questo significa che R(1) e RL occupano la stessa dislocazione dellaparte reale di C mentre R(2) ed IM occupano la parte immaginaria diC. Le istruzioni COMMON ed EQUIVALENCE sono di difficile com-prensione e diventano veramente utili in sistemi composti da moltesubroutines tipo pacchetti grafici o matematici.

1.51 Istruzione EXTERNAL

L’ istruzione EXTERNAL permette che nomi simbolici possano essereusati come argomenti in istruzioni CALL e in chiamate a funzioni, senzache il compilatore crei automaticamente una variabile nel momentodella referenza. In pratica specifica che un nome e un simbolo globaledefinito fuori dell’ unita del programma. Si presentano due casi:

• a) identificare sottoprogrammi o punti ENTRY passati come ar-gomenti,

26

CAPITOLO 1. IL FORTRAN-77 1.52. ISTRUZIONE INTRINSIC

• b) identificare un blocco DATA nel momento del link del pro-gramma.

1.52 Istruzione INTRINSIC

Questa istruzione indica che un certo nome simbolico e una funzioneintrinseca e puo essere passata come argomento in istruzioni CALLe chiamate a funzioni. Il nome della funzione intrinseca specificatomantiene il tipo di variabile associato. Riportiamo un esempio sull’ usodi INTRINSIC/EXTERNAL:♣ INTRINSIC SINH,COSH♣ nna7 EXTERNAL MYFUNC♣ CALL GRAPH (SINH)♣ CALL GRAPH (MYFUNC)♣ END♣ SUBROUTINE GRAPH (FUNC)♣ EXTERNAL FUNC♣ CALL PLOT (FUNC,0.0, 3.14)♣ END

1.53 Istruzione SAVE

Il FORTRAN permette l’ allocazione in memoria sia dinamica che sta-tica. Le variabili, i vettori e i blocchi COMMON dichiarati nell’ unitaMAIN, sono allocati staticamente e mantengono sempre la loro defini-zione. Le variabili, i vettori e i blocchi COMMON dichiarati solo insottoprogrammi vengono allocati dinamicamente al momento dell’ ese-cuzione del sottoprogramma. Quando questo esegue l’ istruzione RE-TURN o END, sono deallocati e perdono la loro definizione(e chiara-mente il loro valore!). L’ istruzione SAVE permette di mantenere la de-finizione di queste variabili e di salvare il loro contenuto, specificandolein due diversi modi:♣ SAVE nome1,nome2♣ SAVE

Nel primo sono salvati nome1 e nome2, mentre nel secondo tuttigli argomenti della sottunita. IN UNIX-FORTRAN questo avviene giaanche senza l’ istruzione SAVE.

27

1.54. ISTRUZIONE DATA CAPITOLO 1. IL FORTRAN-77

1.54 Istruzione DATA

Questa istruzione e adoperata per definire i valori iniziali delle variabilie dei vettori. Strettamente parlando non e un’ istruzione di specificazio-ne, segue tutte le istruzioni di dichiarazione e prende effetto solamentequando il programma e inserito nella memoria. La forma e del tipo:♣ DATA nlista/clista/nlista/clista/

dove ogni nlista e una lista di variabili, vettori, elementi di vettori,sottostringhe oppure liste di DO implicite; ogni clista e una lista dicostanti che forniscono il valore iniziale. Le virgole separanti il simbolo/ e nlista sono degli optional.

I nomi di argomenti fittizi o funzioni non possono essere adoperatiin codesta istruzione. Riportiamo alcuni esempi:♣ PARAMETER (NX=64,NY=80,NTOTAL=NX*NY)♣ REAL*4 ALLSET (NX,NY)♣ DATA ALLSET /NTOTAL * 0.0 /♣ CHARACTER * 8 A,B♣ DATA A,B,C/’VOGLIAMO’/’TRONCARE’/

1.55 Istruzioni di assegnazione

Le istruzioni di assegnazione sono usate per assegnare un valore, usual-mente quello di un’ espressione, ad una variabile o ad un elemento divettore. Possono anche essere utilizzate per modificare i contenuti dilocazioni di memoria assoluta.

1.56 Assegnazioni aritmetiche

Le istruzioni di assegnazione aritmetica sono usate per assegnare unvalore a variabili aritmetiche. Assumono la forma:

v = espr

dove v e una variabile o elemento di vettore di un certo tipo aritmeticoed espr e un’ espressione aritmetica. Se il tipo di espr e diverso daquello di v, allora il valore di espr viene convertito al tipo di v primadella memorizzazione: questo potrebbe causare dei troncamenti.

28

CAPITOLO 1. IL FORTRAN-77 1.57. ASSEGNAZIONI LOGICHE

1.57 Assegnazioni logiche

Le istruzioni di assegnazione logica assegnano un valore ad una variabilelogica ed operano nello stesso modo delle assegnazioni aritmetiche:

v = espr

dove v e una variabile o elemento di vettore di tipo logico ed espr eun’ espressione logica o anche di tipo aritmetico. Questa ultima puocontenere espressioni relazionali o espressioni di tipo CHARACTER.Se il tipo di espr non e logico, il valore assegnato a v e. false., se ilvalore di espr e zero, mentre e. true. se il valore di espr e diverso dazero. Esempio:♣ LOGICAL LIMITE (1:10)♣ REAL*4 R(1:100)♣ LIMITE(J)=R(J).GE.RMIN.AND.R(J).LE.RMAX

1.58 Assegnazioni CHARACTER

Ha la forma del tipo:v = espr

dove v e una variabile o elemento di vettore o ancora sottostringa ditipo CHARACTER ed espr e un’ espressione CHARACTER. Se lalunghezza di v e maggiore di espr vengono aggiunti degli spazi alladestra di espr fino a raggiungere la lunghezza di v; se v e minore di espri caratteri in eccesso sulla destra sono rimossi.

1.59 Istruzione ASSIGN

Questa istruzione e usata per assegnare l’ indirizzo di un istruzioneetichettata ad una variabile intera. La forma e del tipo:♣ ASSIGN S TO V

dove S e un’ etichetta di un istruzione eseguibile o un’ istruzioneFORMAT presente nella stessa unita e V e una variabile intera. Lavariabile V puo essere poi adoperata nello stesso programma sia inun’ istruzione GO TO sia come identificatore di FORMAT in un’ istru-zione READ/ WRITE /PRINT. La variabile non puo ovviamente es-sere richiamata in altre espressioni aritmetiche a meno che non siariassegnata con un valore intero. Esempio:

29

1.60. ISTRUZIONI DI CONTROLLO CAPITOLO 1. IL FORTRAN-77

♣ DIMENSION ARRAY (1:10)♣101 FORMAT (10F10.4)♣ ASSIGN 101 TO MYFORM♣ WRITE ( LP, MYFORM) ARRAY

1.60 Istruzioni di controllo

Le istruzioni di controllo vengono utilizzate, in FORTRAN, per dirige-re il flusso dell’ elaborazione: sono inclusi in questo tipo di istruzionecostrutti per il loop, diramazioni condizionate e incondizionate, strut-ture decisionali a scelta multipla e possibilita di pausa o interruzionedell’ esecuzione di un programma. I trasferimenti di controllo in unaunita del programma si effettuano principalmente con l’ uso di istru-zioni GO TO e IF e, inoltre, con l’ uso sporadico di END ed ERR inistruzioni di input/output.

1.61 Istruzione GO TO incondizionata

Questa istruzione produce un trasferimento incondizionato del flusso suun’ altra istruzione etichettata:♣ GO TO S

dove S e un etichetta di un’ istruzione eseguibile presente nella stessaunita.

1.62 Istruzione GO TO calcolata

Questa istruzione trasferisce il controllo ad una istruzione scelta da unalista secondo una condizione specifica. La forma e del tipo:♣ GO TO(s1,s2,....sn) espr

dove ogni sn e un etichetta di un istruzione eseguibile presente nellastessa unita ed espr e un’ istruzione di tipo intero utilizzata per selezio-nare,tra le sn, l’ istruzione cui trasferire il controllo. La virgola primadi espr e un optional. La scelta viene effettuata nel modo seguente: seil valore di espr e 1, il controllo viene trasferito ad s1, se espr ha valore2, il controllo viene trasferito ad s2, e cosı via. Se il valore di espr eminore di 1 o maggiore del numero di etichette nella lista, non viene

30

CAPITOLO 1. IL FORTRAN-77 1.63. ISTRUZIONE GO TO ASSEGNATA

effettuato alcun trasferimento. La stessa etichetta puo comparire piudi una volta nella stessa lista. Esempio:♣ GO TO(100,200,300,...100) MOD(KCONTA,5)+1

1.63 Istruzione GO TO assegnata

E un’ istruzione obsoleta ma mantenuta per compatibilita con il FOR-TRAN IV. Trasferisce il controllo ad una istruzione che ha un’ etichettaassegnata ad un intero tramite un ASSIGN. Ha la forma:♣ GO TO V

dove V e la variabile intera cui e stato assegnato l’ indirizzo di unaistruzione eseguibile etichettata.

1.64 Blocco IF

Un blocco IF permette di eseguire le istruzioni in maniera condizionale.La forma e del tipo:♣ IF (LEXP1) THEN♣ BLOCCO1♣ ELSEIF (LEXP2) THEN♣ BLOCCO2♣ ELSE♣ BLOCCO3♣ END IF

dove ogni LEXP e un’ espressione logica e ogni blocco contiene istru-zioni eseguibili (eccetto l’ END), compresi altri blocchi IF. Ogni bloccodi istruzioni che segue un IF / ELSEIF e eseguito solo se la corrispon-dente espressione logica ha il valore .true.; il blocco ELSE e eseguitosolo se tutte le precedenti espressioni logiche sono risultate false. Dopol’ esecuzione delle istruzioni contenute in un blocco, il controllo vienetrasferito all’ istruzione immediatamente seguente l’ ENDIF. Le istru-zioni ELSE ed ELSE IF possono essere omesse: l’ ultima in particolarepuo essere ripetuta quante volte si vuole.

Per ogni istruzione IF-THEN deve essere sempre presente l’ ENDIF. Ricordiamo che non e possibile passare dal blocco in un’ altra partedell’ unita.

31

1.65. ISTRUZIONE IF LOGICA CAPITOLO 1. IL FORTRAN-77

1.65 Istruzione IF logica

L’ istruzione logica IF permette che una singola istruzione sia eseguibilein maniera condizionale. L’ istruzione assume la forma:♣ IF (LEXP) ISTRUZ

dove LEXP e un’ espressione logica ed ISTRUZ e un istruzioneeseguibile,escluso DO, IF-THEN, ELSE IF, ELSE, END IF, END o unaltro IF logico. L’ istruzione ISTRUZ e eseguita solo se l’ espressionelogica LEXP ha il valore.TRUE.. Esempio:♣ IF (X.GT.0.5 ) CALL MODIFICA (X)

1.66 Istruzione IF aritmetica

Quest’ istruzione e obsoleta e viene mantenuta solo per compatibilitacon ilFORTRAN-IV. La sua forma e del tipo:♣ IF (ESPR)s1,s2,s3

dove ESPR e un’ espressione INTEGER, REAL o DOUBLE PRE-CISION, e ogni s e un etichetta relativa ad una istruzione eseguibile.Il flusso e trasferito ad s1, s2, s3 a seconda che il valore di ESPR sianegativo, zero o positivo.

1.67 Istruzione DO

L’ istruzione DO e usata per definire un DO-LOOP. Questo loop puocontenere un numero a piacere di istruzioni eseguibili ( eccetto l’ END )che vengono eseguite ripetutamente con una sequenza regolare di valoridi una variabile chiamata variabile del do. La forma dell’ istruzione edel tipo:♣ DO s, v=e1,e2,e3

oppure♣ DO s, v=e1,e2, e

dove s e l’ etichetta di un’ istruzione eseguibile, v e la variabiledel DO di tipo INTEGER, REAL o DOUBLE PRECISION e ogni ee un’ espressione di tipo intero, reale o doppia precisione, indicanti,rispettivamente, il valore iniziale, il valore finale e l’ incremento di v.La virgola dopo l’ etichetta e un optional ma e raccomandato controgli errori. Esempio:

32

CAPITOLO 1. IL FORTRAN-77 1.68. ISTRUZIONE DO-WHILE

♣ DO 55, N= -100, MAX (LIMIT, 100),2♣ blocco istruzioni♣55 continue

Il DO-LOOP consiste nel blocco di istruzioni eseguibili dall’ istru-zione DO a quella con etichetta s inclusa. Il LOOP puo contenereistruzioni CALL, funzioni di referenza e blocchi IF/ END IF come al-tri blocchi DO. E necessario pero che ogni struttura, sia blocco IF cheDO-LOOP, venga interamente contenuta nel range del DO-LOOP piuesterno. In alcune implementazioni, tipo UNIX-FORTRAN o MICRO-SOFT, e stata soppressa l’ etichetta del DO e abbiamo semplicemente:♣ DO v=e1, e2, e3♣ blocco istruzioni♣ END DO

L’ esecuzione di un DO-LOOP procede nel modo seguente:

1. viene valutato e1 (valore iniziale) ed assegnato alla variabile v delloop, con appropriata conversione di tipo, se necessaria.

2. vengono valutate le espressioni e2 ed e3 ( rispettivamente valo-re limite ed incremento o passo ). Se e3 e stato omesso, vieneconsiderato il valore di default 1.

3. Viene calcolato il numero di iterazioni con la seguente espressione:MAX ( INT (( e2 - e1 + e3 )/e3 ), 0 ). Si ottiene quindi il numerodi volte che il blocco di istruzioni nel loop verra eseguito.

4. se questo numero di istruzioni e nullo, il controllo dell’ esecuzio-ne viene trasferito alla prima istruzione che segue l’ END DO ol’ etichetta del loop.

5. se il numero di iterazioni e maggiore di zero, vengono eseguite leistruzioni del blocco.

6. la variabile del loop viene incrementata di e3, il numero di itera-zioni e decrementato di uno e si riprende dal passo 4.

1.68 Istruzione DO-WHILE

Questa e un istruzione che esiste solo sul UNIX-FORTRAN o MICRO-SOFT e fornisce un metodo per il loop non necessariamente controllato

33

1.69. ISTRUZIONE CONTINUE CAPITOLO 1. IL FORTRAN-77

da un conto di iterazioni. In pratica si esegue un blocco di istruzionifino a che il valore dell’ espressione logica e falsa. In ogni passo delLOOP viene controllata la variabile logica; se e vera (.true. ) vengonoeseguite le espressioni del LOOP, se e falsa (.false. ) il controllo vienetrasferito all’ istruzione dopo il LOOP. La forma dell’ istruzione e laseguente:♣ [ DO s ] WHILE expr

dove s e l’ istruzione etichettata eseguibile che definisce il ran-ge del loop: puo essere anche una istruzione REPEAT o END DO.L’ espressione expr e un test di tipo logico. Il seguente esempio e trattodalla teoria dei frattali e costituisce il loop principale dell’ insieme diMandelbrot:♣ IM = 0♣ DO I= 0,NPIX1♣ CY=YMIN + I*DELTA/NPIX♣ DO J=0,NPIX1♣ CX=XMIN+J*DELTA/NPIX♣ COMPLEX = CX + IMM*CY♣ C DIMENSIONE NUMERO COMPLESSO♣ Z=0♣ NDIM = 0♣ DO WHILE (CABS(Z).LE. 2.0.AND. NDIM .LT.ITER)♣ Z=Z*Z+COMPLEX♣ NDIM=NDIM+1♣ END DO♣ IM = IM +1♣ MATRIX (IM) = NDIM♣ END DO♣ END DO

1.69 Istruzione CONTINUE

L’ istruzione CONTINUE fornisce un punto di riferimento e, pratica-mente, traferisce il controllo all’ istruzione seguente. E usato princi-palmente come fine del vecchio DO LOOP, ma puo comparire ovunquenell’ unita di programma. Esempi:♣ GO TO 400♣ ...........

34

CAPITOLO 1. IL FORTRAN-77 1.70. ISTRUZIONE STOP

♣400 CONTINUE♣ DO 600 J=1,300♣ ..................♣600 CONTINUE

1.70 Istruzione STOP

Questa istruzione termina l’ esecuzione del programma e produce unritorno al sistema operativo usato. Si presentano due casi:♣ STOP

oppure♣ STOP [disp]

dove [disp] rappresenta una variabile CHARACTER oppure un nu-mero non piu grande di 5 cifre. Questa variabile puo essere un messaggioda visualizzare o un numero di riferimento per errori di esecuzione delprogramma. Esempio:♣ STOP ’Ci sono troppi elementi’

1.71 Istruzione PAUSE

L’ istruzione PAUSE e simile all’ istruzione dello STOP con la differenzache la sospensione dell’ esecuzione e solo temporanea. L’ esecuzioneriprende con una istruzione del sistema operativo che nel caso del VMSe $ CONTINUE. Si sconsiglia l’ uso di questa istruzione quando i sistemidevono essere trasportabili. La forma e del tipo♣ PAUSE [s]

dove s e una stringa di caratteri o un numero fino a 5 cifre, inviatasu uno schermo come messaggio per l’ utente.

1.72 Istruzioni di INPUT/OUTPUT

Le istruzioni di Input/Output forniscono un canale di comunicazionetra il programma FORTRAN e il mondo esterno. Ci sono tre gruppidi istruzioni di input/output. Nel primo gruppo si trovano le istruzioniREAD, WRITE, PRINT che attuano il trasferimento dei dati fra unita;fanno parte del secondo gruppo BACKSPACE, ENDFILE e REWIND

35

1.73. NOZIONI GENERALI CAPITOLO 1. IL FORTRAN-77

che compiono operazioni ausiliarie. Tutte queste istruzioni fanno riferi-mento ad unita logiche intere che si riferiscono a file ed apparati esterni( lettori di nastri / stampanti / dischetti / terminali ). Nel terzo gruppotroviamo OPEN, CLOSE e INQUIRE che correlano le unita logiche diuscita o ingresso con il nome dei file.

1.73 Nozioni generali

In questo paragrafo verranno spiegati in generale i concetti che riguar-dano record e file per una migliore comprensione delle istruzioni diInput / Output.

1.74 Record

Tutti gli I/O del FORTRAN hanno luogo tramite una struttura chiama-ta RECORD. Un record puo essere un singolo carattere, una sequenzadi caratteri o di valori; puo essere una linea di testo, le coordinate peril pennino di un plotter, i dati letti da uno scanner o ( fino a circa 10anni fa ) una scheda perforata.

Il FORTRAN usa tre tipi di record: formattato, non formattato e difine file (” endfile ” ). Un record formattato e costituito da una sequenzadi caratteri ASCII, terminata o meno con un ” carriage return ”, unavanzamento di linea o entrambi. Ogni singola riga di testo di questomanuale, per esempio, e un record formattato. La lunghezza minimaper un record formattato e zero e la massima e 1024. Un record nonformattato e una sequenza di valori e la sua interpretazione dipende daltipo di dati. Per esempio il binario 01010111 puo essere interpretatocome valore intero 87 o valore carattere ” W ” a seconda del tipo didato. La lunghezza minima di un record non formattato e zero e lamassima e 1024, eccetto per i record ad accesso sequenziale ( vedi oltre) non contenenti informazioni sulla lunghezza del record, che puo quindiessere illimitata. Il record di endfile e l’ ultimo record di un file e nonha lunghezza.

36

CAPITOLO 1. IL FORTRAN-77 1.75. FILE

1.75 File

Un file e composto da record, puo essere creato ed e accessibile anchecon mezzi diversi dai programmi in linguaggio FORTRAN. Per esempiosi puo creare e modificare un file di testo con un editor e manipolarne leinformazioni con un programma in FORTRAN. I file contenuti in nastrie dischi sono generalmenti detti esterni; i file memorizzati nella memoriaprincipale sono chiamati file interni. La posizione all’ interno di un filesi riferisce al prossimo record che sara letto o scritto. Quando si accedead un file ( lo si ” apre ” ) si e posizionati prima del primo record e lafine del file e subito dopo il record di endfile. Alcune istruzioni di I/Opermettono di cambiare la posizione corrente in un file.

1.76 Accesso ai file

Il metodo usato per trasferire record a o da file e chiamato modalita diaccesso. L’ accesso ai file avviene attraverso il loro nome. I file esternipossono contenere sia record formattati che non formattati. Quando unrecord in un file puo essere letto o scritto in modo arbitrario, casuale,la modalita di accesso e detta DIRETTA. Si accede ai singoli recordattraverso un numero di record, che e un intero positivo. Tutti i recordin un file ad accesso diretto hanno la stessa lunghezza e contengono soloi dati attualmente scritti su di essi; non ci sono caratteri di fine record.I record possono essere riscritti ma non cancellati. Generalmente soloi file su disco possono usare la modalita di accesso diretto per il trasfe-rimento dei record. Quando i record vengono trasferiti in ordine, unodopo l’ altro, la modalita di accesso e detto SEQUENZIALE. I recordin un file ad accesso sequenziale possono essere di diversa lunghezza.Alcuni file, come terminali, stampanti e lettori di nastri, possono usaresolo questa modalita di accesso. I file ad accesso sequenziale formattatocontengono solitamente informazioni di tipo testo e ogni record ha uncarattere di fine record. L‘ accesso sequenziale non formattato e gene-ralmente utilizzato per due motivi che possono sembrare contrastanti:1) per il controllo di plotters, terminali grafici, stampanti, etc. e perprocedere su informazioni binarie non codificate come i file.OBJ. Inquesto caso e importante che i dati trasferiti da e al mezzo esterno sia-no un flusso di bytes senza informazioni sulla lunghezza del record. 2)per la comprensione dei dati e la velocita di accesso che caratterizzano

37

1.77. FILE INTERNI CAPITOLO 1. IL FORTRAN-77

questa modalita. In questo caso deve poter essere possibile determinarela lunghezza del record per letture parziali.

1.77 File interni

I file interni possono contenere variabili CHARACTER, elementi divettore di tipo CHARACTER e sottostringhe CHARACTER. Un fi-le interno che contenga una variabile, un elemento di vettore o unasottostringa, tutti di tipo CHARACTER, ha un solo record, la cui lun-ghezza e quella dell’ entita CHARACTER considerata. Un file internoche contenga un vettore di tipo CHARACTER, ha tanti record quantisono gli elementi del vettore, la cui lunghezza e uguale a quella deglielementi stessi. I dati di questi file possono essere trasferiti solo conla modalita di accesso sequenziale formattato. Solitamente questi filevengono utilizzati per la conversione del tipo di variabili da numericoa CHARACTER.

1.78 Istruzione OPEN

Questa istruzione e adoperata per connettere un’ unita logica ad un filee specificare le sue caratteristiche. Puo aprire un file esistente oppurecrearne uno nuovo. Assume la forma del tipo♣ OPEN ( olista )

dove OLISTA e una lista di argomenti ognuno separato da unavirgola:

UNIT= u Specifica il numero dell’ unita. I numeri 5, 6, (*)indicano lo standard input ( tastiera ) e lo standard output (schermo ).

STATUS = cexpr Dove cexpr e un’ espressione CHARACTERche assume il valore ’OLD’ se il file e gia esistente, ’NEW’ se voglia-mo crearne uno nuovo, ’SCRATCH’ per inserire un file nella memoriadinamica.

FILE = cexpr dove cexpr e una espressione CHARACTER cheriporta il nome del file; non deve essere usata insieme a ’SCRATCH’.

ACCESS = cexpr dove cexpr e un’ espressione CHARACTERche assume il valore ’SEQUENTIAL’ o ( per default ) ’DIRECT’.

FORM = cexpr dove cexpr e un’ espressione CHARACTER conil valore ’FORMATTED’ oppure ’UNFORMATTED’

38

CAPITOLO 1. IL FORTRAN-77 1.79. ISTRUZIONE CLOSE

RECL = iexp dove iexp e un’ espressione intera che specifica lalunghezza del record; deve essere data per i file ad accesso DIRECT.Le unita sono CHARACTERS per i file FORMATTED.

BLANK = cexpr , dove cexpr e un’ espressione CHARACTERche specifica il trattamento standard degli spazi: puo essere un ’NULL’ (gli spazi vengono ignorati ) oppure uno ’ZERO’ ( gli spazi sono sostituitida zeri ).

ERR = s dove s e una etichetta che si riferisce ad una istruzionecui viene passato il controllo dell’ esecuzione in caso di errore.

IOSTAT=ios dove ios e il nome simbolico di una variabile IN-TEGER*4 o di un elemento di vettore. Nel caso di errori, ios assumeun valore positivo; se si verifica la condizione di endfile, ios assume ilvalore -1 ed infine in caso di operazione conclusa senza errori ios assumeil valore zero. Esempi:♣ OPEN (unit=8, file=’espe.dat’, status=’NEW’)♣ OPEN (unit=N+5, file=dev//dati, status=’OLD’,♣ $ ERR = 99 )♣ OPEN (unit= 1,file=fnome, status=’OLD’,♣ $ FORM=’UNFORMATTED’, IOSTAT= KODE )

1.79 Istruzione CLOSE

Tutti i file e le unita sono chiuse automaticamente alla fine del pro-gramma. Una istruzione CLOSE esplicita puo servire per riassegnaresuccessivamente un file o una unita;♣ CLOSE ( lista )

dove lista e una serie di argomenti scelti da:

UNIT = u specifica l’ unita

STATUS = cexp e un’ espressione CHARACTER del tipo ’KEEP’per conservare i file oppure ’DELETE’ per cancellarli.

IOSTAT = iv ritorna il tipo di errore.

ERR = s se si verifica un errore di chiusura, ritorna all’ unita s.

Esempio:♣ CLOSE (UNIT=15, STATUS=’DELETE’,ERR=’999’)

39

1.80. ISTRUZIONE INQUIRE CAPITOLO 1. IL FORTRAN-77

1.80 Istruzione INQUIRE

L’ istruzione INQUIRE e usata per ottenere informazioni sulle proprietadi unita logiche o di file. Ci sono due forme:♣ INQUIRE (unit = u, ilist )♣ INQUIRE (file= cexp, ilist )

dove ilist e una lista contenente uno o piu dei seguenti argomenti:ERR = s specifica un label per errore di uscita NUMBER= iv ri-torna il numero dell’ unita logica RECL = iv ritorna la lunghezza delrecord se il file e ad accesso diretto IOSTAT= iv ritorna lo stato delcodice NEXTREC =iv ritorna il numero del prossimo record da leg-gere o da scrivere EXIST = lv ritorna.true. oppure.false., se l’ unitao file esiste o meno OPENED = lv ritorna.true. oppure.false., sel’ unita o file e connesso o meno NAMED = lv ritorna.true. se ilfile ha un nome, oppure.false. NAME = variabile CHARACTERdefinisce il nome del file ACCESS = cv ritorna ’SEQUENTIAL’ op-pure ’DIRECT’ secondo il tipo di accesso definito in OPEN FORM= cv ritorna ’FORMATTED’ oppure ’UNFORMATTED’ secondo iltipo di file connesso BLANK = cv ritorna un ’NULL’ oppure un’BLANK’ SEQUENTIAL =cv ritorna ’YES’ oppure ’NO’oppureindefinito secondo il tipo di accesso al file DIRECT =cv ritorna’YES’ oppure ’NO’oppure indefinito FORMATTED =cv ritorna’YES’ oppure ’NO’ oppure indefinito secondo il tipo di FORM indi-cato UNFORMATTED =cv ritorna ’YES’ oppure ’NO’ oppureindefinito

1.81 Istruzioni READ/WRITE

L’ istruzione READ trasferisce i dati da un file interno od esterno aduna lista di espressioni o elementi di vettori. L’ istruzione WRITEvaluta una lista di espressioni o elementi di vettori e trasferisce i valorisu file interni od esterni. Le istruzioni hanno la forma seguente:♣ READ (unit= 1,FMT = fmt, REC=iexpp,♣ $ ERR =s,END = s,IOSTAT =iv)rlista♣ WRITE (unit= 1,FMT = fmt,REC=iexpp,♣ $ ERR =s,END = s,IOSTAT =iv)wlista

dove unit = u indica l’ unita cui e connesso il file. FMT =fmt specifica il formato in cui leggere o scrivere i dati. REC = rec

40

CAPITOLO 1. IL FORTRAN-77 1.82. ISTRUZIONE REWIND

specifica il record in caso di connessioni con accesso diretto. ERR= s per ridirigere l’ esecuzione in caso di errore nell’ operazione ditrasferimento dati. END = s per ridirigere l’ esecuzione in caso siverifichi la condizione di endfile durante l’ operazione di trasferimentodati. IOSTAT = ios per monitorare gli errori o lo stato di fine file.

Ricordiamo che:

1. il trasferimento di dati ” unformatted ” e permesso solo con fileesterni.

2. viene trasferito un record per ogni istruzione di I/O.

3. il trasferimento di dati ” formatted ” richiede la specifica delformato in cui sono stati scritti o devono essere scritti i dati.

4. per ogni istruzione di I/O possono essere trasferiti uno o piurecord.

L’ istruzione PRINT trasferisce i dati dalla dislocazione interna alterminale in modo interattivo oppure al batch log in modo batch (ossiail comando viene posto in coda ). Prende la forma:♣ Print fmt,list

L’ istruzione TYPE e praticamente uguale al PRINT ma funzio-na solo per alcune implementazioni tipo UNIX-FORTRAN. Assume laforma:♣ TYPE *,list

L’ istruzione ACCEPT trasferisce i dati dal terminale alla memo-ria interna; l’ accesso e di tipo sequenziale. Anche questa istruzionefunziona solo in UNIX-FORTRAN. Riportiamo la forma piu semplice:♣ ACCEPT *,list

Esempi:♣ WRITE (91,101)jx,jy,index(1),index(2)♣ WRITE (5,’(A)’)TEXT

1.82 Istruzione REWIND

Quest’ istruzione e usata per posizionare un file sequenziale all’ iniziodel primo record. Non ha effetto se il file e gia ’riavvolto’. Assume laforma:♣ REWIND ( UNIT = u, ERR =s,IOSTAT =iv)

41

1.83. ISTRUZIONE BACKSPACE CAPITOLO 1. IL FORTRAN-77

Esempio:♣ REWIND ( UNIT = 3,ERR =120 )

1.83 Istruzione BACKSPACE

L’ istruzione puo essere usata su un file sequenziale per muovere laposizione del file indietro di un record. Assume la forma del tipo:♣ BACKSPACE u

dove u e una variabile intera che specifica l’ unita logica interessata,ad esempio un file su disco.

1.84 Istruzione ENDFILE

L’ istruzione ENDFILE e usata per scrivere un record di end alla finedi un file. Questo record in UNIX-FORTRAN corrisponde ad un bytecon il valore ASCII 26 ( CTRL/Z). Prende la forma:♣ ENDFILE u

dove u e una variabile intera che specifica l’ unita logica interessata.Supponiamo adesso di volere aggiungere alcuni dati su un file gia esi-stente. Dovremo prima arrivare alla fine del file tramite un istruzionedi READ, ritornare indietro di una posizione tramite BACKSPACE,inserire tramite WRITE i dati interessati, richiudere il file tramite EN-DFILE e finalmente chiudere il file. Come esempio pratico riportiamola SUBROUTINE SCRIVIDATI:

SUBROUTINE scrividati (ave,sdev )

implicit real *4 (a-h,o-z)

LOGICAL EXST

iunit = 12

iunit2 = 13

OPEN (iunit2,file=’spessore.dat ’,type=’old’)

read (iunit2 , 101) spessore_par,icasi

100 format (f10.4,f10.4,f10.4)

101 format (f10.4,i2)

INQUIRE(FILE=’risultati.dat’,EXIST=EXST )

IF (EXST) THEN

OPEN(iunit ,file=’risultati.dat’,type=’old’)

120 read (iunit , 100,end=130 )zz, xx ,yy

42

CAPITOLO 1. IL FORTRAN-77 1.85. LA SPECIFICAZIONE FORMAT

go to 120

130 continue

backspace iunit

write (iunit , 100 )spessore_par,ave,sdev

endfile iunit

close (iunit )

close (iunit2)

ELSE

OPEN(iunit,file=’risultati.dat’,type=’new’)

write (iunit , 100)spessore_par, ave, sdev

close (iunit )

close (iunit2)

END IF

return

END

1.85 La specificazione FORMAT

La specificazione del formato dei dati e richiesta dalle istruzioni READ/WRITE oppure PRINT per controllarne la traduzione fra forma internaed esterna. Un FORMAT puo essere inserito in qualsiasi parte delprogramma ed e caratterizzato da una etichetta relativa posta nelleprime cinque colonne. Assume la forma:♣100 FORMAT (ed1,ed2.....)

dove ed1, ed2...edn sono i cosidetti descrittori di formato che sidividono come segue.

1.86 Descrittori ripetibili

Ogni descrittore di questo tipo corrisponde ad una operazione di input/output. Riportiamo i piu usati:

IW In input, trasferisce ”w” caratteri dal campo esterno ai va-lori di lista e assegna ad essi un valore intero. Il dato esterno deveessere un intero. In output trasferisce il corrispondente valore ad uncarattere esterno che e lungo ”w” elementi. Fw.d In input trasferisce”w” caratteri dal campo esterno e assegna ad essi, come valore reale, icorrispondenti elementi di lista ( che devono essere reali ). Se il primo

43

1.86. DESCRITTORI RIPETIBILI CAPITOLO 1. IL FORTRAN-77

Tabella 1.3: Tabella sul format Gw.d

Magnitudine del dato Conversione effettiva

m < 0.1 Ew.d[Ee]0.1 ≤ m < 1.0 F (w − 4).d, n(′ ′)1.0 ≤ m < 10.0 F (w − 4).(d− 1), n( ′ ′)

. .

. .

. .10 ∗ ∗d− 2 ≤ m < 10 ∗ ∗d− 1 F (w − 4).1, n( ′ ′)

10 ∗ ∗d− 1 ≤ m < 10 ∗ ∗d F (w − 4).0, n( ′ ′)m ≥ 10 ∗ ∗d Ew.d[Ee]

carattere non nullo e un segno -, il campo e trattato come negativo. Seinvece il campo non contiene ne un numero decimale ne un esponenteviene trattato come un numero reale di ”w” cifre in cui le ”d” cifre sonoa destra rappresentate da degli zeri. In output, trasferisce la corrispon-dente lista arrotondata alla ”d” posizione decimale. Il termine ”w”deve essere grande abbastanza per includere il segno ( se necessario ),almeno una cifra alla sinistra del punto decimale e ”d” cifre alla destradel decimale. Ew.d In input, funziona come il formato F. In output,E trasferisce il valore della corrispondente lista, arrotondata a ”d” cifredecimali con ”E” che specifica il valore dell’ esponente. Il termine wdeve essere grande abbastanza per includere: segno,uno zero,un puntodecimale, ”d” cifre e un esponente. Quindi per tenere conto di questeesigenze il termine ”w” deve essere maggiore o uguale a d+ 7.

Dw.d Vale per i reali in doppia precisione. In output si comportacome il formato E, eccetto che la lettera D rimpiazza la lettera E prece-dente l’ esponente e abbiamo due cifre come esponente. Gw.d In input,funziona come il formato F. In output trasferisce la lista corrisponden-te in funzione della magnitudine, vedi tabella seguente. Tenendo contodelle varie possibilita ”w” deve essere piu grande o uguale a ”d”+8.

nHc1..cn In input trasferisce ”n” caratteri dal campo esterno al de-scrittore medesimo.Il primo carattere appare subito dopo la H. Prendeil nome di variabile di Hollerith. Lw In input trasferisce ”w” caratteridal campo esterno e assegna un valore logico al corrispondente elemen-

44

CAPITOLO 1. IL FORTRAN-77 1.87. DESCRITTORI NON RIPETIBILI

to di lista. Se il primo carattere non nullo e T,t,.T,.t viene assegnatoil valore.TRUE. al corrispondente elemento se invece i primo caratterenon nullo vale F, f,.F, or.f, viene assegnato il valore.FALSE.. In outputtrasferisce la lettera T ( se la variabile e vera ) o la lettera F ( se lavariabile e falsa ). La lettera T o F viene messa nella posizione piua destra preceduta da w-1 spazi. Aw Se il corrispondente elementodi I/O e di tipo CHARACTER viene trasmesso un dato character. Ilvalore di ”w” deve essere minore od uguale a 32767. Il ”w” puo essereomesso e per default viene assegnato un numero di caratteri uguali aquelli dichiarati nella variabile numerica.

1.87 Descrittori non ripetibili

I seguenti descrittori servono a controllare la posizione in un record o inun file, ma non corrispondono ad una operazione di I/O: ’stringa’ dovela stringa viene mandata in uscita Tn si posiziona tramite il TAB allacolonna n Trn si posiziona tramite il TAB a n colonne a destra Tln siposiziona tramite il TAB a n colonne a sinistra nX lascia n spazi SPproduce un segno prima di tutti i numeri seguenti positivi SS sopprimeil segno prima di tutti i numeri seguenti positivi BN tratta i blanksseguenti in input come nulli BZ tratta i blanks seguenti in input comezeri KP setta il fattore di scala per i susseguenti I/O / comincia unnuovo record

Facciamo alcuni esempi:♣ A = 5♣ B = - 10♣ C = 2224♣ WRITE ( 5, 20 ) A, B, C♣20 FORMAT(3I4)♣ WRITE ( 5, 21 ) A, B, C♣21 FORMAT(3(1X,I4))♣ WRITE ( 5, 22 ) A, B, C♣22 FORMAT(1X,’A = ’,I2 ,1X,’B = ’,I3,1X,’C = ’,I5,/)♣ A = 12♣ B = 25.4♣ C = -1.23♣ WRITE ( 5, 20 ) A, B, C♣20 FORMAT(F3.0,F4.1,F5.2)

45

1.87. DESCRITTORI NON RIPETIBILI CAPITOLO 1. IL FORTRAN-77

♣ WRITE ( 5, 21 ) A, B, C♣21 FORMAT(F3.0,1X,2(E12.4,1X))♣ WRITE ( 5, 22 ) A, B, C♣22 FORMAT(F3.0,1X,E12.4,1X,D12.4)♣ A = 1.♣ B = 10.♣ C = 100.♣ D = 1000.♣ WRITE ( 5, 22 ) A, B, C, D♣22 FORMAT(5G10.4)♣ WRITE ( 5, 22 ).TRUE.,.FALSE.♣22 FORMAT(2L2)♣ WRITE ( 5, 22 ) ’BUONGIORNO’♣22 FORMAT(A6)

46

Capitolo 2

Il Laboratorio

Nei prossimi capitoli verranno illustrati alcuni programmi per l’ ela-borazione dei dati ricavati dalle esperienze di laboratorio. Questi pro-grammi sono stati sviluppati in ambiente MS-DOS usando il compi-latore Microsoft-FORTRAN versione 5.0. Le esercitazioni sono stateconcepite ed implementate in modo tale da risultare accessibili sin dallaprima volta anche agli utenti sprovvisti di ogni nozione di informatica.L’ utente e infatti guidato attraverso tutte le fasi dell’ esercitazione,dall’ immissione corretta dei dati all’ elaborazione dei parametri sta-tistici di interesse ed infine alla presentazione finale dei risultati daistruzioni che appaiono di volta in volta sullo schermo video o grafi-co. Per alcuni programmi viene dapprima illustrato l’ argomento o gliargomenti a cui esso si riferisce e le formule utilizzate. Viene inoltreinserito il listato del programma in linguaggio FORTRAN,commentatonelle varie fasi, per l’ utente che desidera apportare modifiche parti-colari. I programmi sono quasi tutti scritti in ambiente MS-DOS inFORTRAN-77 e fanno uso di subroutines contenute nel libro [REC 92]e per semplicicita non contengono uscite grafiche.

2.1 RETTA

Questo programma illustra come, date due serie di datix ey e l’ er-rore relativo alla variabile y in interattivo,al fine di stabilire eventualirelazioni di tipo lineare tra di esse, vengono calcolati:

47

2.1. RETTA CAPITOLO 2. IL LABORATORIO

1. valor medio:

x =1

N

N∑i=1

xi (2.1)

y =1

N

N∑i=1

yi (2.2)

2. deviazione standard:

sx =

√√√√ 1

N − 1

N∑i=1

(xi − x)2 (2.3)

sy =

√√√√ 1

N − 1

N∑i=1

(yi − y)2 (2.4)

3. coefficienti della retta di regressione (coef. angolare b e intercettaa):

b =

∑Ni=1(xi − x)(yi − y)∑N

i=1(xi − x)2(2.5)

a = y − bx (2.6)

4. errore sui coefficienti

error(b) =

√√√√ ∑Ni=1(yi − a− bxi)2

(N − 2)∑Ni=1(xi − x)2

(2.7)

error(a) =

√√√√∑Ni=1(yi − a− bxi)2

(N − 2)[

1

N+

x2∑Ni=1(xi − x)2

] (2.8)

5. coefficiente di correlazione r di Pearson

r =

∑Ni=1(xi − x)(yi − y)√∑N

i=1(xi − x)2∑Ni=1(yi − y)2

(2.9)

6. Valore di χ2

χ2(a, b) =N∑i=1

((yi − a− bxi)

σi)2 (2.10)

48

CAPITOLO 2. IL LABORATORIO 2.1. RETTA

Segue il listato del programma che usa le subroutine FIT e PEARSNestratte dal libro [REC 92]:

program retta

C definizione delle variabili e dei vettori utilizzati nel programma

implicit real (a-h,o-z)

character ris

parameter ( max =100)

dimension xreal (1:max)

dimension yreal (1:max)

dimension sigma (1:max)

print *,’Programma per il calcolo statistico ’

print *,’Media,Deviazione standard e correlazione lineare ’

print *,’Calcolo del chiquadro e del coefficiente ’

print *,’di correlazione lineare ’

call attesa

print *,’riportiamo anche il parametro q che indica la bonta’

print *,’del fit , vedi Numerical Recipes I/II ’

print *,’Q = gammq ( (N-2)/2 , chi**2 /2 ) ’

print *,’N ==> numero di dati ’

print *,’dove gammq e la funzione gamma incompleta ’

print *,’ Q > 0.1 ==> fit accettabile ’

print *,’ Q > 0.001 ==> fit accettabile solo se gli ’

print *,’ errori sono non-normali oppure ’

print *,’ sottostimati ( vedi modulo Youg ’

print *,’Q < 0.001 ==> il modello non funziona !!! ’

call attesa

100 format (a)

print *,’Hai messo i dati nel file retta.dat a ’

print *,’triplette x,y errore sulle y ? ’

print *,’rispondere Y/N ? ’

read 100,ris

if ((ris.eq.’Y’).or.(ris.eq.’y’)) then

go to 150

elseif ((ris.eq.’N’).or.(ris.eq.’n’)) then

stop ’Metti i dati nel file retta.dat a triplette x,y,sy !!’

else

stop ’SCEMO , risposta non valida !!!!!’

end if

49

2.2. PARABO CAPITOLO 2. IL LABORATORIO

150 continue

call attesa

open (unit=90,file=’retta.dat’, status=’old’)

C input dati

n = 0

199 continue

read (90,*,end = 200)xx,yy ,sy

n = n + 1

xreal ( n ) = xx

yreal ( n ) = yy

sigma ( n ) = sy

go to 199

200 continue

close ( unit = 90 )

mwt = 1

call fit (xreal,yreal,n,sigma,mwt,a,b,errora,errorb,chi2,q)

call pearsn (xreal,yreal,n,r,prob,z)

print *,’fit della retta y = a + b * x ’

print *,’a =’,a

print *,’b =’,b

print *,’errore su a = ’,errora

print *,’errore su b = ’,errorb

print *,’valore del chiquadro =’,chi2

print *,’valore del coefficiente lineare r =’,r

print *,’valore di bonta del fit Q =’,q

stop

end

2.2 PARABO

Questo programma viene utilizzato per l’ analisi dei dati dell’ esercita-zione sulla bilancia di Cavendish (determinazione della costante di gra-vitazione mediante il metodo della parabola) oppure per l’ esercitazionesulla rotaia ad aria (determinazione della accelerazione di gravita). Sidevono immettere il numero n di dati e, successivamente le n triplettedi valori x,y e errore su y. Tramite la subroutine LFIT, ricavata dallibro [REC 92] , si calcolano i coefficienti a,b,c della curva parabolica di

50

CAPITOLO 2. IL LABORATORIO 2.2. PARABO

regressione y = a + b× x + c × x2 e i relativi errori con il metodo deiminimi quadrati. Viene anche visualizzato il valore sperimentale dellavariabile χ2. Segue il listato del programma.

program parabo

C

C definizione delle variabili e dei vettori

C utilizzati nel programma

C

implicit real (a-h,o-z)

parameter ( max =50 )

parameter ( nterm = 3 )

dimension x (1:max)

dimension y (1:max)

dimension sig (1:max)

dimension a (1:nterm)

dimension covar (1:nterm,1:nterm)

dimension lista (1:nterm)

external funcs

C

print *,’Programma per il calcolo statistico ’

print *,’Fit con un polinomio di secondo grado’

99 continue

C

C input dei dati da tastiera

C

call attesa

print *,’Avete messo i dati nel file parabo.dat ?’

print *,’In sequenza , x , y , errore sulle y .....’

call attesa

open (unit=11,file=’parabo.dat’,status=’old’)

n = 0

199 continue

read (11,*,end = 200)xx,yy,errory

n = n + 1

x ( n ) = xx

y ( n ) = yy

sig ( n ) = errory

go to 199

51

2.3. PENDO CAPITOLO 2. IL LABORATORIO

200 continue

close ( unit = 11)

do 201 jj = 1, nterm

lista ( jj ) = jj

201 continue

call lfit (x,y,sig,n,a,lista,nterm,covar,nterm,chisq,funcs)

errora = sqrt ( covar (1,1))

errorb = sqrt ( covar (2,2))

errorc = sqrt ( covar (3,3))

print *,’fit del polinomio y = a + b * x + c * x**2 ’

print *,’a =’,a(1)

print *,’b =’,a(2)

print *,’c =’,a(3)

print *,’errore su a = ’,errora

print *,’errore su b = ’,errorb

print *,’errore su c = ’,errorc

print *,’Valore del chiquadro =’,chisq

stop

end

2.3 PENDO

Questo programma e utile nell’ analisi dei dati dell’ esercitazione con ilpendolo composto. Vengono richiesti gli errori sulle misure della varia-bile y (tempo) per la prima e la seconda serie di dati. Successivamentisi immettono il numero di dati e le coppie di valori x (distanza),y (tem-po) per la prima serie di dati. Il programma valuta immediatamentei coefficienti della curva parabolica y = a + b × x + c × x2, median-te chiamata alla subroutine LFIT estratta dal libro [REC 92] che usail metodo dei minimi quadrati quando il polinomio interpolante e disecondo grado.

Analogo procedimento viene seguito per la seconda serie di dati. Ilprogramma calcola poi l’intersezione delle parabole, l’ errore sull’inter-sezione e finalmente il valore di g ed il suo errore. Alleghiamo il listatodel programma.

program pendo

implicit real (a-h,o-z)

parameter ( max =50 )

52

CAPITOLO 2. IL LABORATORIO 2.3. PENDO

parameter ( nterm = 3 )

parameter ( inc = 3 )

parameter ( reg = 980.5343)

character * 1 prova,video,rispo

character * 40 messag

dimension x (1:max)

dimension y (1:max)

dimension sig (1:max)

dimension x2 (1:max)

dimension y2 (1:max)

dimension sig2 (1:max)

dimension result(1:20)

dimension a (1:nterm)

dimension covar (1:nterm,1:nterm)

dimension lista (1:nterm)

dimension para (1: inc )

dimension epara (1: inc )

dimension messag (1: inc )

external funcs

prova =’y’

el = 99.3

errl = 0.1

teta = 5.0

eteta = 1.0

nosci = 30

if (( prova.eq.’Y’).or.(prova.eq.’y’)) then

print *,’Programma per il calcolo di g in CGS tramite ’

print *,’intersezione di due parabole ’

print *,’ovviamente deve essere presente l’’andamento parabolico’

print *,’nei dati ..... ’

print *,’i dati x rappresentano le distanze della massa mobile’

print *,’dal coltello e sono espresse in cm ’

print *,’i dati y rappresentano i tempi impiegati a compiere ’

print *,’n oscillazioni a vanno espressi in secondi ’

print *,’gli errori sulle y (errori tempi) pure vanno espressi ’

print *,’in secondi ! ’

99 continue

100 format (a1)

print *,’Hai messo i dati nel file pendo1.dat ’

53

2.3. PENDO CAPITOLO 2. IL LABORATORIO

print *,’a triplette x [cm],y [sec] errore sulle y [sec] ? ’

print *,’rispondere Y/N ? ’

read 100, rispo

if ((rispo.eq.’Y’).or.(rispo.eq.’y’)) then

go to 150

elseif ((rispo.eq.’N’).or.(rispo.eq.’n’)) then

stop ’Metti i dati nel file pendo1.dat a triplette x,y,sy !!’

else

stop ’SCEMO , risposta non valida !!!!!’

end if

150 continue

print *,’Hai messo i dati nel file pendo2.dat ’

print *,’a triplette x [cm],y [sec] errore sulle y [sec] ?’

print *,’rispondere Y/N ! ’

read 100, rispo

if ((rispo.eq.’Y’).or.(rispo.eq.’y’)) then

go to 151

elseif ((rispo.eq.’N’).or.(rispo.eq.’n’)) then

stop ’Metti i dati nel file pendo2.dat a triplette x,y,sy !!’

else

stop ’SCEMO , risposta non valida !!!!!’

end if

151 continue

print*,’Inserire l(distanza fra i coltelli)in [cm] tipo 99.3’

read *,el

print*,’Inserire errore l in [cm] tipo 0.1 ’

read *,errl

print*,’Inserire il numero di oscillazioni tipo 30 ’

read *,nosci

print*,’Inserire ampiezza iniziale [gradi] tipo 5 ’

read *,teta

print*,’Inserire errore ampiezza iniziale [gradi] tipo 1 ’

read *,eteta

call attesa

end if

open (unit=90,file=’pendo1.dat’, status=’old’)

open (unit=91,file=’pendo2.dat’, status=’old’)

n = 0

199 continue

54

CAPITOLO 2. IL LABORATORIO 2.3. PENDO

read (90,*,end = 200)xx,yy ,sy

n = n + 1

x ( n ) = xx

y ( n ) = yy

sig ( n ) = sy

go to 199

200 continue

n2 = 0

299 continue

read (91,*,end = 300)xx,yy ,sy

n2 = n2 + 1

x2 ( n2 ) = xx

y2 ( n2 ) = yy

sig2 ( n2 ) = sy

go to 299

300 continue

errort = sig2 (1)

close ( unit = 90 )

close ( unit = 91 )

mfit = nterm

do 333 j=1,mfit

lista(j)= j

333 continue

call lfit (x,y,sig,n,a,lista,nterm,covar,nterm,chisq,funcs)

errora = sqrt ( covar (1,1))

errorb = sqrt ( covar (2,2))

errorc = sqrt ( covar (3,3))

print *,’Risultati prima parabola ....................... ’

print *,’Minimi quadrati con y = a + b * x + c * x**2 ’

print *,’numero dati =’,n

print *,’a =’,a(1)

print *,’b =’,a(2)

print *,’c =’,a(3)

print *,’errore su a = ’,errora

print *,’errore su b = ’,errorb

print *,’errore su c = ’,errorc

print *,’Valore del chiquadro =’,chisq

if (( prova.eq.’Y’).or.(prova.eq.’y’)) then

call attesa

55

2.3. PENDO CAPITOLO 2. IL LABORATORIO

end if

result (1) = a(1)

result (2) = a(2)

result (3) = errora

result (4) = errorb

result (5) = chisq

result (6) = a(3)

result (7) = errorc

mfit = nterm

do 334 j=1,mfit

lista(j)= j

334 continue

call lfit (x2,y2,sig2,n2,a,lista,nterm,covar,nterm,chisq,funcs)

errora = sqrt ( covar (1,1))

errorb = sqrt ( covar (2,2))

errorc = sqrt ( covar (3,3))

print *,’Risultati seconda parabola ....................... ’

print *,’Minimi quadrati con y = a + b * x + c * x**2 ’

print *,’numero dati =’,n

print *,’a =’,a(1)

print *,’b =’,a(2)

print *,’c =’,a(3)

print *,’errore su a = ’,errora

print *,’errore su b = ’,errorb

print *,’errore su c = ’,errorc

print *,’Valore del chiquadro =’,chisq

if (( prova.eq.’Y’).or.(prova.eq.’y’)) then

call attesa

end if

result (11) = a(1)

result (12) = a(2)

result (13) = errora

result (14) = errorb

result (15) = chisq

result (16) = a(3)

result (17) = errorc

a1 = result (1)

b1 = result (2)

c1 = result (6)

56

CAPITOLO 2. IL LABORATORIO 2.3. PENDO

a2 = result (11)

b2 = result (12)

c2 = result (16)

call radici (a1,b1,c1,a2,b2,c2,val,err,errort)

val = val / nosci

err = err / nosci

para (1) = val

para (2) = el

para (3) = teta

epara (1) = err

epara (2) = errl

epara (3) = eteta

messag (1)= ’peso-errore periodo =’

messag (2)= ’peso-errore lunghezza =’

messag (3)= ’peso-errore ampiezza =’

print *,’Dati di input per il calcolo di g ’

print *,’periodo come da intersezione [sec]=’,para(1)

print *,’errore periodo [sec]=’,epara(1)

print *,’distanza fra i coltelli [cm ]=’,para(2)

print *,’errore distanza fra i coltelli[cm ]=’,epara(2)

print *,’ampiezza iniziale [gradi ]=’,para(3)

print *,’errore ampiezza iniziale [gradi ]=’,epara(3)

print *,’************************************** ’

print *,’Dati di output ’

print*,’g=4*(pigrec**2)*(l/T**2)(1+(0.5*sen(0.5*teta))**2)**2’

call value (valore,para,inc)

error = 100 * abs ((reg - valore)/reg)

vero = valore

print *,’g [cm/sec**2] =’,valore

print *,’g vero a TORINO [cm/sec**2]=’,reg

video=’N’

call errori (valore,evalo,para,inc,epara,messag,video)

ippm = nint ( evalo /vero * 1.0E+6 )

print * ,’errore su g [cm/sec**2] =’,evalo

print * ,’errore su g in ppm =’,ippm

print * ,’errore g sul valore vero =’,error,’ %’

teta=0.

para (3) = teta

call value (valoresenza,para,inc)

57

2.3. PENDO CAPITOLO 2. IL LABORATORIO

ippmteta=nint(((valoresenza-vero)/vero)* 10E+6)

print *,’g senza teta =’,valoresenza

print *,’influenza di teta su g in ppm =’,ippmteta

video=’Y’

valore = vero

call errori (valore,evalo,para,inc,epara,messag,video)

stop

end

subroutine radici (a1,b1,c1,a2,b2,c2,val,err,errort)

implicit real ( a-h,o-z)

c = a1 - a2

b = b1 - b2

a = c1 - c2

disc = b**2 - 4 *a*c

if ( disc.lt.0.0) then

stop’radici immaginarie’

end if

x1 = - ( b/(2*a)) + sqrt ( disc)/(2*a)

x2 = - ( b/(2*a)) - sqrt ( disc)/(2*a)

y1 = a1 + b1 * x1 + c1 *x1**2

y2 = a1 + b1 * x2 + c1 *x2**2

if ((y1.gt.0).and.(y2.gt.0)) then

print *,’caso di due radici positive ’

val = ( y1 + y2 )/2.0

err = abs ( y1-y2)

if ( err.lt.errort) then

err = errort

end if

end if

if ((y1.lt.0).and.(y2.gt.0)) then

print *,’caso di una radice positiva e una negativa’

val = y2

err = errort

end if

if ((y2.lt.0).and.(y1.gt.0)) then

print *,’caso di una radice positiva e una negativa’

val = y1

err = errort

end if

58

CAPITOLO 2. IL LABORATORIO 2.3. PENDO

return

end

subroutine value (valore,para,inc)

implicit real (a-h,o-z )

parameter ( pigrec = 3.1415927 )

dimension para ( 1 : inc )

t = para (1)

el = para (2)

teta = para (3)

tetar = 2 * pigrec * teta / 360.

cor = ( 1.0 + (0.5 * sin (tetar/2.0))**2.0)**2.0

valore = (4 * pigrec**2 * el /( t **2))*cor

return

end

subroutine errori (valore,evalo,para,inc,epara,messag,video)

implicit real (a-h,o-z)

character * 40 messag

character * 1 video

parameter ( int = 10 )

parameter ( norma=10000 )

dimension para (1:inc )

dimension epara (1:inc )

dimension work (1:int )

dimension square (1:int )

dimension inter (1:int )

dimension messag (1:inc )

if (int.lt.inc)stop ’cambiare int’

vero = valore

evalo = 0.0

do 100 j = 1 ,inc

work (j) = para (j)

100 continue

quadro = 0.0

do 200 j = 1 ,inc

do 190 i = 1 , inc

para (i) = work(i)

190 continue

59

2.4. CAVEN CAPITOLO 2. IL LABORATORIO

para (j) = para(j) + epara(j)

call value (valore,para,inc)

delta = valore - vero

d2 = delta *delta

quadro = quadro + d2

square (j) = d2

200 continue

sq2 = 0.0

do 205 j= 1,inc

sq2 = sq2 + square ( j)

205 continue

do 206 j= 1,inc

inter (j) = nint ( square ( j) * (norma*1.0)/sq2)

206 continue

if (( video.eq.’Y’).or.( video.eq.’y’)) then

print*,’************************************’

print *,’Peso normalizzato a’,norma,’ del’’errore delle variabili’

do 207 j= 1, inc

print *,messag(j),inter(j)

207 continue

end if

evalo = sqrt ( quadro)

do 210 i = 1 , inc

para (i) = work(i)

210 continue

return

end

2.4 CAVEN

Questo programma valuta i parametri caratteristici dell’ oscillatoresmorzato, con il quale vengono interpolati i dati delle ampiezze di oscil-lazione [cm] in funzione del tempo [sec] ricavati dalla esperienza conla bilancia di Cavendish. Le coppie di dati (tempi ed ampiezze) vannomemorizzate nel file CAV.DAT, mentre in fase di esecuzione viene ri-chiesto l’ errore [cm] sulla misura delle ampiezze e la distanza [cm] fraasse della bilancia e riga per la rilevazione dati. Tramite la subrouti-ne MRQMIN (estratta dal libro [REC 92]) vengono calcolati i cinque

60

CAPITOLO 2. IL LABORATORIO 2.4. CAVEN

parametri A,T,α, τ ,Cost che caratterizzano la seguente equazione:

y = A× seno(2× π × t/T + α)× e−tτ + Cost (2.11)

In seguito si ricava la costante di Cavendish utilizzando i seguenti valoricaratteristici della bilancia:

m = (15.0± 0.12)× 10−3 Kg (2.12)

M = (1.503± 0.001) Kg (2.13)

r = (4.5± 0.1)× 10−2 m (2.14)

d = (5.0± 0.1)× 10−2 m (2.15)

l ≈ 5.0m (misurato durante l′ esercitazione) (2.16)

Riportiamo la formula utilizzata per il calcolo;

G = A×B × C ×D (2.17)

dove

A =d× r2

2M(2.18)

B = ω2 + τ−2 dove ω =2× πT

(2.19)

C = |ε2 − ε1| (2.20)

con

ε1 =E1

lE1 = posizione iniziale sull′ asta graduata (2.21)

ε2 =E2

lE2 = posizione finale sull′ asta graduata (2.22)

D =1

1− r3

[(2d)2 + r2]32

(2.23)

Viene inoltre calcolato l’ errore percentuale sul valore esatto della co-stante dato da

Errore = 100× |Gricavato −Gvero|Gvero

. (2.24)

61

2.4. CAVEN CAPITOLO 2. IL LABORATORIO

Viene calcolato poi l’ errore sulla costante G usando la legge di pro-pagazione degli errori di Gauss: in questo caso le derivate parzialivengolo calcolate tramite algoritmo numerico. Il programma calcolapoi G tramite il metodo dell’ escursione finale e quello della parabola.Per confronto viene poi riportato Riportiamo il listato del programmaCAVEN:

program caven

implicit real (a-h,o-z )

character * 1 prova,video,rispo

parameter(max=250,ma =5,inc=7 )

parameter(nterm = 2 )

dimension x (1:max )

dimension y (1:max )

dimension sig (1:max )

dimension a (1:ma )

dimension ia (1:ma )

dimension covar(1:ma,1:ma)

dimension alpha(1:ma,1:ma)

dimension erroa (1:ma )

dimension para (1:inc )

dimension epara (1:inc )

dimension ovar (1:nterm,1:nterm)

dimension lista (1:nterm)

external fun2c

external fosci

prova =’N’

1000 format (A2)

100 format (a1)

print *,’Programma per il calcolo statistico ’

print *,’Interpolazione di una sinusoide smorzata ’

print *,’*************************************************’

if ((prova.eq.’N’).or.(prova.eq.’n’)) then

print *,’Avete messo i dati nel file caven.dat ? ’

print *,’in sequenza t [sec] , ampiezze [cm] ? ’

print *, ’inserire (Y/N ) ! ’

read 100, rispo

if ((rispo.eq.’Y’).or.(rispo.eq.’y’)) then

go to 150

62

CAPITOLO 2. IL LABORATORIO 2.4. CAVEN

elseif ((rispo.eq.’N’).or.(rispo.eq.’n’)) then

stop ’Metti i dati nel file caven.dat a coppie !!’

else

stop ’SCEMO , rispo non valida !!!!!’

end if

150 continue

end if

print *,’Attenzione il programma non funziona sempre e quindi’

print *,’osservate le seguenti regole !!!!!!!!!!! ’

print *,’1)il primo dato rappresenta inizio della sinusoide’

print *,’2)il dato ultimo deve essere quello di equilibrio finale’

print *,’3)i dati anomali vanno eliminati ’

call attesa

if ((prova.eq.’Y’).or.(prova.eq.’y’)) then

el = 544.0

eel = 1.0

eey = 0.2

else

print *,’Scrivere il valore della distanza [cm] ( circa 544)’

read *,el

print *,’Scrivere errore della distanza [cm] ( circa 1 )’

read *,eel

print *,’Scrivere errore lettura scala [cm] ( circa 0.2)’

read *,eey

end if

open (unit=90,file=’caven.dat’, status=’old’)

npt = 0

199 continue

read (90,*,end = 200)xx,yy

npt = npt + 1

if ( npt.gt.max) stop ’Aumentare max !’

x ( npt ) = xx

y ( npt ) = yy

sig ( npt ) = eey

go to 199

200 continue

close ( unit = 90 )

a (1) = abs ( y (npt) - y (1))

a (2) = 500

63

2.4. CAVEN CAPITOLO 2. IL LABORATORIO

a (3) = asin (( y (1)- y (npt))/a(1))

a (4) = 500

a (5) = y (npt)

do 210 j= 1 , ma

ia (j) = 1

210 continue

alamda=-1

call mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,

* MA,chisq,fosci,alamda)

k=1

itst=0

1 continue

k=k+1

ochisq=chisq

call mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,

* MA,chisq,fosci,alamda)

if (chisq.gt.ochisq) then

itst=0

else if (abs(ochisq-chisq).lt.0.1) then

itst=itst+1

endif

if (itst.lt.4) then

goto 1

endif

alamda=0

call mrqmin(x,y,sig,NPT,a,ia,MA,covar,alpha,

* MA,chisq,fosci,alamda)

erroa ( 1 ) = sqrt(covar (1,1))

erroa ( 2 ) = sqrt(covar (2,2))

erroa ( 3 ) = sqrt(covar (3,3))

erroa ( 4 ) = sqrt(covar (4,4))

erroa ( 5 ) = sqrt(covar (5,5))

print *,’Fittiamo i dati con la seguente funzione nonlineare’

print *,’y=a(5)+a(1)*(sin(2*pigrec*x/a(2)+a(3)))*exp(-x/a(4))’

print *,’a(1) =’,a(1),’ err a(1)=’,erroa(1)

print *,’a(2) =’,a(2),’ err a(2)=’,erroa(2)

print *,’a(3) =’,a(3),’ err a(3)=’,erroa(3)

print *,’a(4) =’,a(4),’ err a(4)=’,erroa(4)

print *,’a(5) =’,a(5),’ err a(5)=’,erroa(5)

64

CAPITOLO 2. IL LABORATORIO 2.4. CAVEN

print *,’chisq =’,chisq

call attesa

recav = 6.67259 E-11

egran = 1.500

er = 4.652 E-02

eb = 4.500 E-02

ed = 5.0 E-02

eegran = 1.0 E-3

eer = 0.1 E-02

eed = 0.1 E-02

para(1) = egran

para(2) = el

para(3) = er

para(4) = ed

para(5) = a(1)

para(6) = a(2)

para(7) = a(4)

epara(1) = eegran

epara(2) = eel

epara(3) = eer

epara(4) = eed

epara(5) = erroa(1)

epara(6) = erroa(2)

epara(7) = erroa(4)

print *,’Sinopsi dei dati fisici adoperati ’

print *,’M [kg] =’,egran,’+-’,eegran

print *,’r [cm] =’,er,’+-’,eer

print *,’b [cm] =’,eb,’+-’,eer

print *,’d [cm] =’,ed,’+-’,eed

print *,’l [cm] =’,el,’+-’,eel

call attesa

call gtutto (cav,para,inc)

ver = cav

valore = cav

video =’N’

cav = ver

call calco (ecav,cav,para,inc,epara,video)

cav = ver

error = 100 * abs ((recav - cav)/recav)

65

2.4. CAVEN CAPITOLO 2. IL LABORATORIO

ippm = nint ( ecav/cav * 1.0E+6 )

print *,’Metodo delle oscillazioni smorzate completo ’

print *,’Formule come in Guida alla Elaborazione dei dati ’

print *,’ oppure Guida alla Fisica del Laboratorio ’

print * ,’costante di Cavendish =’,valore

print * ,’errore costante di Cavendish =’,ecav

print * ,’errore costante di Cavendish in ppm=’,ippm

print * ,’errore percentuale sul valore vero=’,error,’%’

video =’Y’

cav = ver

call calco (ecav,cav,para,inc,epara,video)

call attesa

para ( 3 ) = eb

call gescur (cave,caves,para,inc)

valore = cave

video =’N’

cave = valore

call calcos (ecave,cave,para,inc,epara,video)

cave = valore

ippm = nint ( ecave/cave * 1.0E+6 )

error2 = 100 * abs ((recav - cave )/recav)

print *,’Metodo escursione finale ( come da manuale Leybold )’

print *,’ pigrec **2 b**2 d S ’

print *,’ G = ------------------------------ ’

print *,’ M1 T**2 L ’

print *,’ ’

print *,’ b** 3 ’

print *,’ beta = -------------------------------------- ’

print *,’ (b**2 + 4 d**2 ) sqrt (b**2 + 4 d**2 ) ’

print *,’Costante di Cavendish senza correzione = ’,caves

print *,’Costante di Cavendish con correzione = ’,valore

print *,’Errore Costante con correzione = ’,ecave

print *,’Errore costante con correzione in ppm =’,ippm

print * ,’errore percentuale sul valore vero=’,error2,’%’

video =’Y’

cave = valore

call calcos (ecave,cave,para,inc,epara,video)

call attesa

tquar = a (2) / 2.5

66

CAPITOLO 2. IL LABORATORIO 2.4. CAVEN

do 500 j= 1, npt

if (x( j ).gt.tquar) then

jvec = j -1

go to 501

end if

500 continue

501 continue

do 600 jj = 1 , nterm

lista ( jj ) = jj

600 continue

C LAVORO

print *,’Metodo dell’’ accelerazione ’

print *,’Formule come in Guida alla Fisica del Laboratorio ’

call lfit (x,y,sig,jvec,a,lista,nterm,ovar,nterm,chisq,fun2c)

errora = sqrt ( ovar (1,1))

errorc = sqrt ( ovar (2,2))

print *,’fit della parabola y = a + c * x**2 ’

print *,’tempo massimo =’,tquar

print *,’elementi selezionati =’,jvec

print *,’a =’,a(1)

print *,’c =’,a(2)

print *,’errore su a =’,errora

print *,’errore su c =’,errorc

print *,’Valore del chiquadro =’,chisq

aa = 2 * a (2)

print *,’accelerazione=>2*c[cm/sec**2]=’,aa

call attesa

para (1) = egran

para (2) = el

para (3) = er

para (4) = ed

para (5) = aa

epara(1) = eegran

epara(2) = eel

epara(3) = eer

epara(4) = eed

epara(5) = errorc * 2.0

inc2 = 5

call gpara (cav,para,inc2)

67

2.4. CAVEN CAPITOLO 2. IL LABORATORIO

error = 100 * abs ((recav - cav)/recav)

ippm = nint ( ecav/cav * 1.0E+6 )

vcav = cav

video =’N’

cav = vcav

call capa2 (ecav,cav,para,inc2,epara,video)

cav = vcav

error = 100 * abs ((recav - vcav)/recav)

ippm = nint ( ecav/vcav * 1.0E+6 )

print *,’Metodo della parabola ’

print *,’formule come su Guida alla Fisica del laboratorio ’

print *,’ b**2 d a ’

print *,’ G = -------------- ’

print *,’ 2 M L ’

print *,’ ’

print *,’Costante di Cavendish =’,vcav

print *,’errore costante di Cavendish =’,ecav

print *,’errore costante di Cavendish in ppm=’,ippm

print *,’errore percentuale sul valore vero=’,error,’%’

video =’Y’

cav = vcav

call capa2 (ecav,cav,para,inc2,epara,video)

stop

end

SUBROUTINE fosci(x,a,y,dyda,na)

implicit real (a-h,o-z)

parameter ( pigrec = 3.1415927)

dimension a (1:na)

dimension dyda(1:na)

y=a(5)+a(1)*(sin(2*pigrec*x/a(2)+a(3 )))* exp(-x/a(4))

dyda(1) = (sin(2*pigrec*x/a(2)+a(3 )))* exp(-x/a(4))

ccc = 2*pigrec*x/(a(2)**2.0)

bbb = (cos(2*pigrec*x/a(2)+a(3 )))* exp(-x/a(4))

aaa = - a(1)

dyda(2) = aaa * bbb * ccc

dyda(3) = a(1)*(cos(2*pigrec*x/a(2)+a(3 )))* exp(-x/a(4))

bbb = a(1)*(sin(2*pigrec*x/a(2)+a(3 )))* exp(-x/a(4))

aaa = x /(a(4)**2.0)

dyda(4) = aaa * bbb

68

CAPITOLO 2. IL LABORATORIO 2.4. CAVEN

dyda(5) = 1.0

return

END

subroutine gtutto (cav,para,inc)

implicit real (a-h,o-z)

dimension para (1:inc )

pigrec = 3.1415927

egran = para(1)

el = para(2)

er = para(3)

ed = para(4)

a1 = (ed * er**2)/(2.0* egran )

a2 = (2 * pigrec / para(6))**2 + ( 1/para(7))**2

diff = abs ( atan ( para(5) /el ))

aldif = diff/2.0

a3 = aldif

elle = sqrt ( (2* ed)**2 + er **2 )

a4 = 1.0/(1.0 - (er/elle)**3 )

cav = a1 * a2 * a3 * a4

return

end

subroutine calco (ecav,cav,para,inc,epara,video)

implicit real (a-h,o-z)

character * 1 video

parameter ( int = 10 )

dimension para (1:inc )

dimension epara (1:inc )

dimension work (1:int )

dimension square (1:int )

dimension inter (1:int )

if (int.lt.inc)stop ’cambiare int’

vero = cav

do 100 j = 1 ,inc

work (j) = para (j)

100 continue

quadro = 0.0

do 200 j = 1 ,inc

do 190 i = 1 , inc

para (i) = work(i)

69

2.4. CAVEN CAPITOLO 2. IL LABORATORIO

190 continue

para (j) = para(j) + epara(j)

call gtutto (cavmod,para,inc)

delta = cavmod - vero

quadro = quadro + delta *delta

square ( j )= delta * delta

200 continue

sq2 = 0.0

do 205 j= 1,inc

sq2 = sq2 + square ( j)

205 continue

do 206 j= 1,inc

inter (j) = nint ( square ( j) * 1000./sq2)

206 continue

if (( video.eq.’Y’).or.( video.eq.’y’)) then

print *,’Peso normalizzato a 1000 dell’’errore delle variabili’

print *,’sull’’errore del valore di G-metodo completo ’

print *,’errore su M ( massa grande) =’,inter(1)

print *,’errore su L ( distanza indice) =’,inter(2)

print *,’errore su r ( distanza masse ) =’,inter(3)

print *,’errore su d ( manubrio/2 ) =’,inter(4)

print *,’errore su a(1)(ampiezza oscillazioni) =’,inter(5)

print *,’errore su a(2)(periodo oscillazioni) =’,inter(6)

print *,’errore su a(4)(smorzamento oscillazioni)=’,inter(7)

end if

ecav = sqrt ( quadro)

do 210 i = 1 , inc

para (i) = work(i)

210 continue

return

end

subroutine gescur (cave,caves,para,inc)

implicit real (a-h,o-z)

dimension para (1:inc )

pigrec = 3.1415927

C egran = para(1)

C el = para(2)

C er = para(3)

70

CAPITOLO 2. IL LABORATORIO 2.4. CAVEN

C ed = para(4)

C a(1) = para(5)

C a(2) = para(6)

C a(4) = para(7)

egran = para(1)

el = para(2)

er = para(3)

ed = para(4)

aa1 = pigrec **2

aa2 = er **2

aa3 = ed

aa4 = para(5) * 0.01

dd1 = egran

dd2 = para(6) ** 2

dd3 = el / 100.

ages = aa1 * aa2 * aa3 *aa4 /( dd1 * dd2 * dd3 )

d1 = 1 + 4 * ( ed / er )**2

d2 = sqrt ( 1 + 4 * ( ed / er )**2 )

beta = 1 /( d1 * d2 )

piu = beta * ages

cave = ages + piu

caves = ages

return

end

subroutine calcos (ecave,cave,para,inc,epara,video)

implicit real (a-h,o-z)

character * 1 video

parameter ( int = 10 )

dimension para (1:inc )

dimension epara (1:inc )

dimension work (1:int )

dimension square (1:int )

dimension inter (1:int )

if (int.lt.inc)stop ’cambiare int’

vero = cave

do 100 j = 1 ,inc

work (j) = para (j)

71

2.4. CAVEN CAPITOLO 2. IL LABORATORIO

100 continue

quadro = 0.0

do 200 j = 1 ,inc

do 190 i = 1 , inc

para (i) = work(i)

190 continue

para (j) = para(j) + epara(j)

call gescur (cmod,caves,para,inc)

delta = cmod - vero

quadro = quadro + delta *delta

square (j) = delta * delta

200 continue

sq2 = 0.0

do 205 j= 1,inc

sq2 = sq2 + square ( j)

205 continue

do 206 j= 1,inc

inter (j) = nint ( square ( j) * 1000./sq2)

206 continue

if (( video.eq.’Y’).or.( video.eq.’y’)) then

print *,’Peso normalizzato a 1000 dell’’errore delle variabili’

print *,’sull’’errore del valore di G-escursione finale’

print *,’errore su M ( massa grande) =’,inter(1)

print *,’errore su L ( distanza indice) =’,inter(2)

print *,’errore su r ( distanza masse ) =’,inter(3)

print *,’errore su d ( manubrio/2 ) =’,inter(4)

print *,’errore su a(1)(ampiezza oscillazioni) =’,inter(5)

print *,’errore su a(2)(periodo oscillazioni) =’,inter(6)

print *,’errore su a(4)(smorzamento oscillazioni)=’,inter(7)

end if

ecave = sqrt ( quadro)

do 210 i = 1 , inc

para (i) = work(i)

210 continue

return

end

SUBROUTINE fun2c(X,AFUNC,MA)

DIMENSION AFUNC(MA)

AFUNC(1)= 1.0

72

CAPITOLO 2. IL LABORATORIO 2.4. CAVEN

AFUNC(2)= X**2.0

END

subroutine gpara (cav,para,inc)

implicit real (a-h,o-z)

dimension para (1:inc )

egran = para(1)

el = para(2)

er = para(3)

ed = para(4)

aa = para(5)

cav = ( aa* ed * er**2 )/(2.0 * egran * el )

return

end

subroutine capa2 (ecav,cav,para,inc,epara,video)

C lavoro

implicit real (a-h,o-z)

character * 1 video

parameter ( int = 10 )

dimension para (1:inc )

dimension epara (1:inc )

dimension work (1:int )

dimension square (1:int )

dimension inter (1:int )

if (int.lt.inc)stop ’cambiare int’

vero = cav

do 100 j = 1 ,inc

work (j) = para (j)

100 continue

quadro = 0.0

do 200 j = 1 ,inc

do 190 i = 1 , inc

para (i) = work(i)

190 continue

para (j) = para(j) + epara(j)

call gpara(cavmod,para,inc)

delta = cavmod - vero

quadro = quadro + delta *delta

73

2.5. GAUSS CAPITOLO 2. IL LABORATORIO

square ( j )= delta * delta

200 continue

sq2 = 0.0

do 205 j= 1,inc

sq2 = sq2 + square ( j)

205 continue

do 206 j= 1,inc

inter (j) = nint ( square ( j) * 1000./sq2)

206 continue

if (( video.eq.’Y’).or.( video.eq.’y’)) then

print *,’Peso normalizzato a 1000 dell’’errore delle variabili’

print *,’sull’’errore del valore di G-metodo parabola ’

print *,’errore su M ( massa grande) =’,inter(1)

print *,’errore su L ( distanza indice) =’,inter(2)

print *,’errore su r ( distanza masse ) =’,inter(3)

print *,’errore su d ( manubrio/2 ) =’,inter(4)

print *,’errore su a ( accelerazione ) =’,inter(5)

end if

ecav = sqrt ( quadro)

do 210 i = 1 , inc

para (i) = work(i)

210 continue

return

end

2.5 GAUSS

Questo programma riguarda l’ esercitazione sul viscosimetro ma puoessere usato per altri scopi. Si parte con una serie di dati inseritanel file GAUSS.DAT che vengono lette dal programma e inserite nelvettore DATI. La subroutine ISTOTUTTO provvede poi a calcolare ilmassimo e il minimo e, dato in modo interattivo il numero di classi,la distribuzione in frequenza. Si ottengono poi i seguenti parametristatistici caratterizzanti la distribuzione sperimentale:

1. valor medio dell’ istogramma

x =1

N

N−1∑i=1

(xi)(yi) (2.25)

74

CAPITOLO 2. IL LABORATORIO 2.5. GAUSS

2. deviazione standard dell’ istogramma

s2 =

√√√√ 1

N − 1

N−1∑i=1

yi × (xi − x)2 (2.26)

3. momento di terzo ordine e asimmetria

s3 =

∑N−1i=1 yi × (xi − x)3

(N − 1)× s23(2.27)

4. momento di quarto ordine e curtosi

s4 =

∑N−1i=1 yi × (xi − x)4

(N − 1)× s24− 3 (2.28)

5. χ2:

χ2 =N−1∑i=1

(yi − fi)2

fi(2.29)

Riportiamo il listato del programma GAUSS:

program gauss

C

C definizione delle variabili e dei vettori utilizzati dal programma

C

implicit real (a-h,o-z)

parameter ( max = 1000)

parameter ( pigrec = 3.141592653)

character * 2 risp

character * 1 pausa

dimension xreal (1:max)

dimension iyreal (1:max)

dimension dati (1:max)

open (unit=90,file=’gauss.dat’, status=’old’)

C

C input dei dati da tastiera

C

print *,’Programma per il calcolo statistico ’

75

2.5. GAUSS CAPITOLO 2. IL LABORATORIO

print *,’Analisi di dati con distribuzione Gaussiana ’

print *,’Avete messo i dati nel file gaussiana.dat ? ’

print *, ’se no uscite !!!!!!!!!!!!!!! ’

print *,’<CR> per continuare !’

read 1000,pausa

ndati = 0

199 continue

read (90,*,end = 200)dato

ndati = ndati + 1

dati ( ndati ) = dato

go to 199

200 continue

250 continue

print *,’Scrivere il numero di intervalli desiderato ’

read *,ndiv

print *,’va bene questo numero di intervalli ? <CR>/NO ’

read 1000,risp

if ((risp.eq.’NO’).or.(risp.eq.’no’)) then

go to 200

end if

jvec = 0

999 format (A1)

1000 format (A2)

ndiv1 = ndiv + 1

call isto (dati,ndati,max,ndiv1,xreal,iyreal,delta)

call mom (ndiv1,xreal,iyreal,xmedio,s2,s3,s4,nsumy)

C calcolo del chi-quadro

chi2 = 0.

do 100 j= 1,ndiv

xmezzo = (xreal(j) + xreal (j+1))/ 2.0

yteo = 1./(sqrt (2.0 * pigrec) *s2 )

yteo = yteo * exp (-((( xmezzo - xmedio)/s2 )**2)/2.0 )

yteo = yteo * nsumy * delta

chi2 = chi2 + (( iyreal (j) - yteo)**2) / yteo

100 continue

print *,’xmedio =’,xmedio

print *,’s2 =’,s2

print *,’s3 =’,s3

print *,’s4 =’,s4

76

CAPITOLO 2. IL LABORATORIO 2.5. GAUSS

print *,’nsumy =’,nsumy

print *,’chi2 =’,chi2

print *,’Vuoi cambiare il numero di intervalli ? Y/N’

read 999,pausa

if ((pausa.eq.’Y’).or.(pausa.eq.’y’)) then

go to 250

end if

stop

end

subroutine isto (dati,ndati,max,ndiv1,xreal,iyreal,delta)

implicit real (a-h,o-z)

dimension dati (1:max)

dimension xreal (1:max)

dimension iyreal(1:max)

C azzeramento variabili

do 101 j= 1 ,ndiv1

xreal (j) = 0

iyreal (j) = 0

101 continue

call maxmin (dati,ndati,xmax,xmin )

xmax = xmax + xmax/1000

ndiv = ndiv1 - 1

delta = ( xmax - xmin )/ndiv

do 102 j = 1 , ndiv1

xreal (j) = xmin + ( j-1 ) * delta

102 continue

do 103 j = 1 , ndati

do 104 i= 1 , ndiv

xl = xreal ( i )

xu = xreal ( i + 1 )

if ( (dati (j).ge.xl).and.(dati(j).lt.xu) ) then

iyreal ( i) = iyreal(i) +1

end if

104 continue

103 continue

return

end

subroutine mom (k,x,ifreq,xmean,s2,s3,s4,n)

implicit real (a-h,o-z)

77

2.5. GAUSS CAPITOLO 2. IL LABORATORIO

dimension x (1:k)

dimension ifreq (1:k)

dimension y (1:1000)

C azzeramento variabili

xmean = 0.

s2 = 0.

s3 = 0.

s4 = 0.

n = 0

C calcolo della frequenza totale

do 100 i=1,k-1

n = n + ifreq (i)

100 continue

C calcolo dei punti centrali

do 101 j=1,k-1

y(j) = ( x(j+1) + x (j) )/2.0

101 continue

C calcolo del valor medio

do 102 j=1,k-1

xmean = xmean + y (j) * ifreq (j)

102 continue

xmean = xmean/ n

C calcolo della deviazione standard

do 103 j=1,k-1

s2 = s2 + ((y(j) - xmean )**2) * ifreq (j)

103 continue

s2 = sqrt (s2/(n-1))

C calcolo della asimmetria

do 104 j =1,k-1

s3 = s3 + ((y(j) - xmean )**3) * ifreq (j)

104 continue

s3 = s3 /((n-1)*( s2**3))

C calcolo della curtosi

do 105 j=1,k-1

s4 = s4 + ((y(j) - xmean )**4) * ifreq (j)

105 continue

s4 = ( s4 /((n-1)*( s2**4))) - 3.

return

end

78

CAPITOLO 2. IL LABORATORIO 2.6. TESTT

2.6 TESTT

Questo programma sostituisce le tabelle relative al test di Student. Siinseriscono in maniera interattiva il numero di gradi di liberta n e laprobabilita considerata. Mediante la subroutine BETAI appartenenteal libro [REC 92] viene poi calcolato il valore critico, xp, associato conil taglio inferiore corrispondente alla probabilita della distribuzione diStudent con n gradi di liberta. Riportiamo il listato del programmaTESTT:

program testt

implicit real (a-h,o-z)

logical condiz

character * 1 pausa

parameter ( tmax = 1000 ,ndiv = tmax * 100)

external betai

print *,’Programma per il calcolo del valore critico nella ’

print *,’distribuzione della variabile casuale t di Student ’

print *,’con nu gradi di liberta , dato alfa ’

print *,’<CR> per continuare !’

read 1000,pausa

999 format (A1)

1000 format (A2)

50 continue

print *,’inserire i gradi di liberta ’

read *, n

100 continue

print *,’immettere alfa, livello di confidenza (0.01,0.05,0.1...)’

read *, alfa

condiz = .false.

delta = tmax / ndiv

do 200 jj= 1, ndiv

t = jj * delta

prob= 1 - betai(0.5*n,0.5,n/(n+t**2))

q = ( 1 + prob ) / 2

qmeno = ( 1 - prob ) / 2

if (.not.condiz) then

79

2.7. TESTC2 CAPITOLO 2. IL LABORATORIO

if ( qmeno.lt. (alfa) ) then

tuno = ( jj-1) * delta

condiz =. true.

end if

endif

if ( qmeno.lt. (alfa/2.0)) then

tdue = ( jj-1) * delta

go to 950

end if

200 continue

stop ’ aumentare tmax !!!’

950 continue

print *,’gradi di liberta‘ =’,n

print *,’livello di significativita‘ =’,alfa

print *,’t critico ad una banda =’,tuno

print *,’t critico a due bande =’,tdue

print *,’Vuoi cambiare le variabili ? Y/N’

read 999,pausa

if ((pausa.eq.’Y’).or.(pausa.eq.’y’)) then

go to 50

end if

stop

end

2.7 TESTC2

Questo programma sostituisce le tabelle della variabile χ2 e usa la sub-routine GAMMQ estratta dal libro [REC 92]. Riportiamo il listato delprogramma TESTC2:

program testc2

implicit real (a-h,o-z)

character pausa

parameter ( c2max = 200 ,ndiv = c2max * 100 )

external gammq

13 continue

print *,’Abbiamo due opzioni ....................’

print *,’fissiamo chi2 e troviamo la probabilita‘ (1)’

80

CAPITOLO 2. IL LABORATORIO 2.7. TESTC2

print *,’fissiamo alfa e troviamo chi2 (2)’

read *,numero

if ( numero.eq.1 ) then

go to 11

elseif( numero.eq.2 ) then

go to 12

else

go to 13

end if

11 continue

print *,’Parte per il calcolo di un integrale ’

print *,’fra chi2 e infinito ’

print *,’della della densita‘ di probabilita‘ CHI-QUADRATO ’

print *,’con nu gradi di liberta ’

call attesa

50 continue

print *,’inserire i gradi di liberta ’

read *, n

100 continue

print *,’immettere il chi2 trovato (0.6,0.9,1.5 ,4...)’

read *,chi2

prob = gammq ( 0.5 * n , 0.5 * chi2 )

print *,’gradi di liberta =’,n

print *,’chi2-impostato =’,chi2

print *,’integrale fra chi2 e infinito =’,prob

print *,’ ’

print *,’Vuoi cambiare le variabili ? Y/N’

read 999,pausa

999 format (a1)

if ((pausa.eq.’Y’).or.(pausa.eq.’y’)) then

go to 50

else

stop

end if

12 continue

print *,’Parte per il calcolo del chi2 corrispondente ’

print *,’ad un livello di significativita‘ alfa ’

print *,’della della densita‘ di probabilita‘ CHI-QUADRATO ’

print *,’con nu gradi di liberta ’

81

2.7. TESTC2 CAPITOLO 2. IL LABORATORIO

call attesa

51 continue

print *,’inserire i gradi di liberta‘ ’

read *, n

print *,’immettere il livello di significativa‘ (0.01,0.05,0.1)’

read *,alfa

delta = c2max / ndiv

if ( n .eq.1 ) then

delta = 0.00001

end if

do 995 jj = 1 , ndiv

chi2 = jj * delta

prob = gammq ( 0.5 * n , 0.5 * chi2 )

alfat = 1 - prob

if ( alfat.gt.alfa ) then

c2sin = chi2p

go to 52

end if

chi2p = chi2

995 continue

52 continue

delta = c2max / ndiv

do 996 jj = 1 , ndiv

chi2 = jj * delta

prob = gammq ( 0.5 * n , 0.5 * chi2 )

if ( prob.lt.alfa ) then

chi2d = chi2p

go to 54

end if

chi2p = chi2

996 continue

54 continue

print *,’chi2 a destra =’,chi2d

print *,’chi2 a sinistra =’,c2sin

print *,’gradi di liberta‘ =’,n

print *,’alfa =’,alfa

print *,’Vuoi cambiare le variabili ? Y/N’

read 999,pausa

if ((pausa.eq.’Y’).or.(pausa.eq.’y’)) then

82

CAPITOLO 2. IL LABORATORIO 2.7. TESTC2

go to 51

else

stop

end if

end

83

2.7. TESTC2 CAPITOLO 2. IL LABORATORIO

84

Capitolo 3

Probabilita e statistica

Adesso cerchiamo di capire come sono scritte le varie subroutine cheabbiamo usato in maniera nascosta perche inserite nel link con le libreriedel libro [REC 92] che sono scritte in FORTRAN-77 standard.

3.1 Gamma

La funzione gamma e definita dall’ integrale

Γ(z) =∫ ∞0

tz−1e−tdt . (3.1)

Quando l’ argomento z e un intero, la funzione gamma diventa il fatto-riale dell’ argomento meno uno:

n! = Γ(n+ 1) . (3.2)

Esistono vari metodi per calcolare la funzione gamma e sembra chequello migliore sia quello derivato da Lanczos. Riportiamo adesso laFUNCTION GAMMLN che dato l’ argomento xx ritorna ln (Γ (xx)):

FUNCTION gammln(xx)

REAL gammln,xx

INTEGER j

DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)

SAVE cof,stp

DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,

*24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,

85

3.2. FATTORIALE CAPITOLO 3. PROBABILITA E STATISTICA

*-.5395239384953d-5,2.5066282746310005d0/

x=xx

y=x

tmp=x+5.5d0

tmp=(x+0.5d0)*log(tmp)-tmp

ser=1.000000000190015d0

do 11 j=1,6

y=y+1.d0

ser=ser+cof(j)/y

11 continue

gammln=tmp+log(stp*ser/x)

return

END

3.2 Fattoriale

Riportiamo la FUNCTION che calcola il fattoriale di n distinguendofra i vari valori e alla fine chiamando la funzione gamma:

FUNCTION factrl(n)

INTEGER n

REAL factrl

CU USES gammln

INTEGER j,ntop

REAL a(33),gammln

SAVE ntop,a

DATA ntop,a(1)/0,1./

if (n.lt.0) then

pause ’negative factorial in factrl’

else if (n.le.ntop) then

factrl=a(n+1)

else if (n.le.32) then

do 11 j=ntop+1,n

a(j+1)=j*a(j)

11 continue

ntop=n

factrl=a(n+1)

else

86

CAPITOLO 3. PROBABILITA E STATISTICA 3.2. FATTORIALE

factrl=exp(gammln(n+1.))

endif

return

END

Un punto interessante e notare che che il fattoriale e giusto pervalori piccoli di n ma e approssimato se passiamo come vedremo adessoal logaritmo del fattoriale. Per quanto riguarda il coefficiente binomialedobbiamo invece dobbiamo passare ai logaritmi per evitare gli overflow.Ricordiamo che il coefficiente binomiale vale:(

n

k

)=

n!

k!(n− k)!0 ≤ k ≤ n , (3.3)

e possiamo introdurre la FUNCTION BICO che calcola il suddettocoefficiente:

FUNCTION bico(n,k)

INTEGER k,n

REAL bico

CU USES factln

REAL factln

bico=nint(exp(factln(n)-factln(k)-factln(n-k)))

return

END

la quale adopera il logaritmo del fattoriale:

FUNCTION factln(n)

INTEGER n

REAL factln

CU USES gammln

REAL a(100),gammln

SAVE a

DATA a/100*-1./

if (n.lt.0) pause ’negative factorial in factln’

if (n.le.99) then

if (a(n+1).lt.0.) a(n+1)=gammln(n+1.)

factln=a(n+1)

else

87

3.3. FUNZIONE BETA CAPITOLO 3. PROBABILITA E STATISTICA

factln=gammln(n+1.)

endif

return

END

3.3 Funzione beta

Passando dagli interi al continuo arriviamo alla funzione beta che cosıdefiniamo

B(z, w) = B(w, z) =∫ 1

0tz−1(1− t)w−1dt , (3.4)

questa e collegata alla funzione gamma dalla relazione

B(z, w) =Γ(z)Γ(w)

Γ(z + w), (3.5)

che e facilmente implementata tramite la FUNCTION beta

FUNCTION beta(z,w)

REAL beta,w,z

CU USES gammln

REAL gammln

beta=exp(gammln(z)+gammln(w)-gammln(z+w))

return

END

3.4 Gamma incompleta

La funzione gamma incompleta e definita come

P (a, x) ≡ γ(a, x)

Γ(a)≡ 1

Γ(a)

∫ x

0ta−1e−tdt (a > 0) . (3.6)

Essa ha come limiti

P (a, 0) = 0 e P (a,∞) = 1 . (3.7)

88

CAPITOLO 3. PROBABILITA E STATISTICA3.4. GAMMA INCOMPLETA

Il complemento di P (a,x) e anche chiamata in maniera confusa funzionegamma incompleta

Q(a, x) ≡ 1− P (a, x) ≡ Γ(a, x)

Γ(a)≡ 1

Γ(a)

∫ ∞x

ta−1e−tdt (a > 0) ,

(3.8)ed ha i seguenti valori limiti

Q(a, 0) = 1 e Q(a,∞) = 0 . (3.9)

La notazione P(a,x),γ(a, x) e Γ(a, x) sono standard; la notazione Q(a,x) la useremo solo in questo contesto. A questo punto riportiamo leFUNCTION che ritornano P(a,x) eQ (a,x ) usando la extra lettera p e q dopo gamm:

FUNCTION gammp(a,x)

REAL a,gammp,x

CU USES gcf,gser

REAL gammcf,gamser,gln

if(x.lt.0..or.a.le.0.)pause ’bad arguments in gammp’

if(x.lt.a+1.)then

call gser(gamser,a,x,gln)

gammp=gamser

else

call gcf(gammcf,a,x,gln)

gammp=1.-gammcf

endif

return

END

FUNCTION gammq(a,x)

REAL a,gammq,x

CU USES gcf,gser

REAL gammcf,gamser,gln

if(x.lt.0..or.a.le.0.)pause ’bad arguments in gammq’

if(x.lt.a+1.)then

call gser(gamser,a,x,gln)

gammq=1.-gamser

else

call gcf(gammcf,a,x,gln)

89

3.4. GAMMA INCOMPLETACAPITOLO 3. PROBABILITA E STATISTICA

gammq=gammcf

endif

return

END

L’ argomento gln e dato tramite serie e procedura delle frazionicontinue:

SUBROUTINE gser(gamser,a,x,gln)

INTEGER ITMAX

REAL a,gamser,gln,x,EPS

PARAMETER (ITMAX=100,EPS=3.e-7)

CU USES gammln

INTEGER n

REAL ap,del,sum,gammln

gln=gammln(a)

if(x.le.0.)then

if(x.lt.0.)pause ’x < 0 in gser’

gamser=0.

return

endif

ap=a

sum=1./a

del=sum

do 11 n=1,ITMAX

ap=ap+1.

del=del*x/ap

sum=sum+del

if(abs(del).lt.abs(sum)*EPS)goto 1

11 continue

pause ’a too large, ITMAX too small in gser’

1 gamser=sum*exp(-x+a*log(x)-gln)

return

END

SUBROUTINE gcf(gammcf,a,x,gln)

INTEGER ITMAX

REAL a,gammcf,gln,x,EPS,FPMIN

PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)

90

CAPITOLO 3. PROBABILITA E STATISTICA 3.5. FUNZIONE ERRORI

CU USES gammln

INTEGER i

REAL an,b,c,d,del,h,gammln

gln=gammln(a)

b=x+1.-a

c=1./FPMIN

d=1./b

h=d

do 11 i=1,ITMAX

an=-i*(i-a)

b=b+2.

d=an*d+b

if(abs(d).lt.FPMIN)d=FPMIN

c=b+an/c

if(abs(c).lt.FPMIN)c=FPMIN

d=1./d

del=d*c

h=h*del

if(abs(del-1.).lt.EPS)goto 1

11 continue

pause ’a too large, ITMAX too small in gcf’

1 gammcf=exp(-x+a*log(x)-gln)*h

return

END

3.5 Funzione errori

La funzione degli errori e la funzione degli errori complementari sonocasi speciali della funzione gamma incompleta e sono quindi ottenutedalle procedure summenzionate. Sono cosı definite

erf(x) =2√π

∫ x

0e−t

2dt , (3.10)

e

erfc(x) ≡ 1− erf(x) =2√π

∫ ∞x

e−t2dt . (3.11)

Esse sono legate alle funzioni gamma incomplete dalle relazioni:

erf(x) = P (1

2, x2) (x ≥ 0) , (3.12)

91

3.6. DISTRIBUZIONI CAPITOLO 3. PROBABILITA E STATISTICA

e

erfc(x) = Q(1

2, x2) (x ≥ 0) . (3.13)

Quindi otteniamo facilmente le FUNCTION ERF e ERFC

FUNCTION erf(x)

REAL erf,x

CU USES gammp

REAL gammp

if(x.lt.0.)then

erf=-gammp(.5,x**2)

else

erf=gammp(.5,x**2)

endif

return

END

FUNCTION erfc(x)

REAL erfc,x

CU USES gammp,gammq

REAL gammp,gammq

if(x.lt.0.)then

erfc=1.+gammp(.5,x**2)

else

erfc=gammq(.5,x**2)

endif

return

END

3.6 Distribuzioni

Analizziamo ora alcune delle distribuzioni di probabilita piu note.

3.6.1 Distribuzione del chiquadro

P (χ2|ν) e definita come la probabilita che il il chi-quadro osservato inun modello corretto sia minore di un certo valore χ2. Il suo comple-mento Q (χ2|ν) e la probabilita che il chi-quadro osservato ecceda ilvalore χ2 tramite chance anche per un modello corretto. In ambedue

92

CAPITOLO 3. PROBABILITA E STATISTICA 3.6. DISTRIBUZIONI

i casi ν e un intero, il numero di gradi di liberta del sistema. Questefunzioni sono cosı collegate alle funzioni gamma incomplete:

P (χ2|ν) = P(ν

2,χ2

2

)= gammp (

ν

2,χ2

2

)(3.14)

Q(χ2|ν) = Q(ν

2,χ2

2

)= gammq (

ν

2,χ2

2

). (3.15)

3.6.2 Distribuzione Beta incompleta

La funzione beta incompleta e definita come

Ix(a, b) ≡Bx(a, b)

B(a, b)≡ 1

B(a, b)

∫ x

0ta−1(1− t)b−1dt (a, b > 0) . (3.16)

Questa funzione ha come valori limite:

I0(a, b) = 0 I1(a, b) = 1 , (3.17)

e vale la relazione di simmetria

Ix(a, b) = 1− I1−x(b, a) . (3.18)

La funzione beta incompleta viene implementata attraverso la FUNC-TION BETAI:

FUNCTION betai(a,b,x)

REAL betai,a,b,x

CU USES betacf,gammln

REAL bt,betacf,gammln

if(x.lt.0..or.x.gt.1.)pause ’bad argument x in betai’

if(x.eq.0..or.x.eq.1.)then

bt=0.

else

bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.-x))

endif

if(x.lt.(a+1.)/(a+b+2.))then

betai=bt*betacf(a,b,x)/a

return

else

betai=1.-bt*betacf(b,a,1.-x)/b

return

endif

END

93

3.6. DISTRIBUZIONI CAPITOLO 3. PROBABILITA E STATISTICA

la quale richiede una routine per valutare la frazione continua:

FUNCTION betacf(a,b,x)

INTEGER MAXIT

REAL betacf,a,b,x,EPS,FPMIN

PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30)

INTEGER m,m2

REAL aa,c,d,del,h,qab,qam,qap

qab=a+b

qap=a+1.

qam=a-1.

c=1.

d=1.-qab*x/qap

if(abs(d).lt.FPMIN)d=FPMIN

d=1./d

h=d

do 11 m=1,MAXIT

m2=2*m

aa=m*(b-m)*x/((qam+m2)*(a+m2))

d=1.+aa*d

if(abs(d).lt.FPMIN)d=FPMIN

c=1.+aa/c

if(abs(c).lt.FPMIN)c=FPMIN

d=1./d

h=h*d*c

aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2))

d=1.+aa*d

if(abs(d).lt.FPMIN)d=FPMIN

c=1.+aa/c

if(abs(c).lt.FPMIN)c=FPMIN

d=1./d

del=d*c

h=h*del

if(abs(del-1.).lt.EPS)goto 1

11 continue

pause ’a or b too big, or MAXIT too small in betacf’

1 betacf=h

return

END

94

CAPITOLO 3. PROBABILITA E STATISTICA 3.7. NUMERI RANDOM

3.6.3 Distribuzione di Student

La distribuzione di Student, denotata con A(t|ν), e utile in parecchicontesti di probabilita, ed in particolare quando si tratta di decidere sedue campioni hanno la stessa media. A(t|ν) e la probabilita, per ν gradidi liberta, che una certa statistica t (misurando la differenza osservatafra medie) sia minore del valore osservato nell’ ipotesi che le mediesiano uguali.Due medie sono significativamente differenti se ad esempioA(t|ν) > 0.99. In altre parole, 1- A(t|ν) e il livello di significativita peril quale l’ ipotesi che le medie siano uguali e scartata. La definizionematematica della funzione e

A(t|ν) =1

ν1/2B(12, ν2)

∫ t

−t

(1 +

x2

ν

)−1 + ν

2 dx . (3.19)

I valori limite sono

A(0|ν) = 0 A(∞|ν) = 1 . (3.20)

A(t|ν) e legata alla funzione beta incompleta Ix (a,b) dalla relazione:

A(t|ν) = 1− I νν+t2

(ν2,1

2

), (3.21)

e quindi adottiamo la funzione betai che abbiamo prima definito.

3.7 Numeri random

Un generatore uniforme e un algoritmo capace di generare dei numerisenza una legge apparente fra 0 ed 1 mentre se vogliamo ottenere nu-meri estratti ad esempio da una distribuzione gaussiana caratterizzatada una certa media e distribuzione standard dobbiamo partire da ungeneratore uniforme e poi modificarlo opportunamente.

3.7.1 Generatori di sistema

Quasi tutti i linguaggi di programmazione contengono una subroutineche con l’ indimenticabile nome ”RAN” produce una serie di numerirandom:

x = ran (iseed) , (3.22)

95

3.7. NUMERI RANDOM CAPITOLO 3. PROBABILITA E STATISTICA

All’ inizio si fissa come primo passo un iseed molto grande e poi siprocede alle chiamate seguenti senza piu modificare iseed. Cambiandoiseed si cambia la sequenza mentre chiamate successive con lo stessoiseed producono sequenze uguali.

Bisogna pero fare attenzione a questi generatori di sistema perchequasi sempre sono dei generatori lineari congruenti che generano dellesequenze di interi I1, I2, I3, ognuno fra 0 ed m-1 ( un numero grande )attraverso la relazione ricorsiva

Ij+1 = aIj + c (mod m) , (3.23)

dove m e chiamato modulo, a moltiplicatore e c incremento. Questa ri-correnza si ripetera eventualmente con un periodo che non e piu grandedi m. Se m, a, c sono scelti in maniera appropriata allora il periodosara massimo ovverosia di lunghezza m. In questo caso tutti gli interifra 0 ed m-1 saranno scelti. La sequenza parte esattamente da questopunto. Il numero reale di ritorno compreso fra 0 ed 1 e generalmenteIj+1/m , cosicche e generalmente minore di 1 ( eccetto che una voltain m chiamate ). Il metodo lineare congruente ha il vantaggio di esseremolto veloce richiedendo solamente poche righe per chiamata e quindie di uso universale. Ha lo svantaggio di non essere libero da correla-zioni sequenziali nelle chiamate successive. Se k numeri random sonochiamati a plottare dei punti in uno spazio k-dimensionale ( con ognicoordinata compresa fra 0 ed 1 ), allora i punti non tendono a coprireuniformemente lo spazio ma piuttosto a distribuirsi su (k-1) piani. Cisaranno almeno circa m1/k piani e se le costanti m, a e c non sono scelteesattamente saranno meno di quelli. Il numero m e usualmente vicino alpiu largo intero rappresentabile ≈ 232. Cosicche ad esempio il numerodi piani sui quali triple di punti giacciono in uno spazio tridimensionalee usualmente non piu grande della radice cubica di 232, circa 1600.

3.7.2 Generatore Random - portatile

Tutto sommato possiamo adoperare il semplice algoritmo moltiplicativo

Ij+1 = aIj (mod m) , (3.24)

se il moltiplicatore a e il modulo m sono scelti con attenzione. I signoriPark e Miller propongono

a = 75 = 16807 m = 231 − 1 . (3.25)

96

CAPITOLO 3. PROBABILITA E STATISTICA 3.7. NUMERI RANDOM

Non e possibile implementare direttamente questo algoritmo poicheil prodotto di a per m eccede il massimo intero possibile e quindidobbiamo passare attraverso una fattorizzazione approssimata di m:

m = aq + r , q = [m/a] , r = m mod a , (3.26)

con le parentesi quadrate indicante la parte intera. Se r e piccolo ed inparticolare r < q, e 0 < z < m-1 puo essere dimostrato che sia a ( zmod q ) e r [ z/ q] stanno nel range 0,..., m-1 e che

az mod m ={a(z mod q)− r[z/q] se e ≥ 0;a(z mod q)− r[z/q] +m altrimenti.

(3.27)

L’ applicazione di questo algoritmo alle costanti summenzionate usai valori q =127773 e r=2836. Questo algoritmo viene implementatoattraverso la FUNCTION RAN0

FUNCTION ran0(idum)

INTEGER idum,IA,IM,IQ,IR,MASK

REAL ran0,AM

PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,

*MASK=123459876)

INTEGER k

idum=ieor(idum,MASK)

k=idum/IQ

idum=IA*(idum-k*IQ)-IR*k

if (idum.lt.0) idum=idum+IM

ran0=AM*idum

idum=ieor(idum,MASK)

return

END

Il periodo di ran0 e 231 -1 ≈ 2.0109. Una peculiarita e che il valore 0non e mai raggiunto al primo tentativo.

3.7.3 Generatore di interi

Il problema consiste nel generare una distribuzione di numeri interi conun processo random fra un numero IMIN e IMAX, frontiere comprese.Il problema e di facile soluzione e consiste nel moltiplicare il numeroreale random compreso fra 0 ed 1 per ∆ = IMAX - IMIN +1, dove

97

3.7. NUMERI RANDOM CAPITOLO 3. PROBABILITA E STATISTICA

l’ addizione di uno significa che vogliamo anche gli interi sulla frontiera.L’ implementazione avviene tramite la FUNCTION G05DYFZAN cheprende il nome dal tentativo di simulare una routine delle NAG.

function g05dyfzan (ilower,iupper,idum)

integer ilower ,iupper,iscelto,idum,g05dyfzan

external ran1

real ran1

delta = iupper - ilower + 1.0

iscelto = int ( ran0 (idum) * delta )

g05dyfzan = ilower + iscelto

if (g05dyfzan .lt. ilower ) stop ’errore g05dyfzan piccolo’

if (g05dyfzan .gt. iupper ) stop ’errore g05dyfzan grande ’

return

end

3.7.4 Distribuzione esponenziale

Nella sezione precedente, abbiamo imparato a generare dei numeri ran-dom con una distribuzione di probabilita uniforme, cosicche la proba-bilita di generare un numero compreso fra x e x + dx, denotata conp(x)dx, e data da

p(x)dx ={dx, e 0 < x<1;0, altrimenti.

(3.28)

La distribuzione e ovviamente normalizzata che significa:∫p(x)dx = 1 . (3.29)

Supponiamo adesso di generare una variabile uniforme x. La distribu-zione di probabilita di y, denotata p(y)dy, e determinata dalla leggefondamentale per le trasformazioni di probabilita, che e semplicemente

|p(y)dy| = |p(x)dx| , (3.30)

oppure

p(y) = p(x)|dxdy| . (3.31)

98

CAPITOLO 3. PROBABILITA E STATISTICA 3.7. NUMERI RANDOM

Assumiamo ora che y (x) ≡ - ln (x), e che p( x) e dato per unadistribuzione uniforme. Avremo allora

p(y)dy = |dxdy|dy = e−ydy , (3.32)

che e distribuita esponenzialmente. Questa distribuzione occorre fre-quentemente nei problemi reali, ad esempio come tempo di attesa fraeventi nei fenomeni governati dalla distribuzione di Poisson. Potete an-che notare che la quantita y

λha la distribuzione di probabilita λe−λy. La

implementazione pratica avviene attraverso la FUNCTION EXPDEV.

FUNCTION expdev(idum)

INTEGER idum

REAL expdev

CU USES ran1

REAL dum,ran1

1 dum=ran1(idum)

if(dum.eq.0.)goto 1

expdev=-log(dum)

return

END

3.7.5 Distribuzione Gaussiana

I metodi di trasformazione possono essere generalizzati a piu di unadimensione. Se x1, x2, x3,.... sono delle variabili random con una pro-babilita di distribuzione congiunta p ( x1, x2, x3,..... )dx1, dx2,dx3...dxne y1,y2,.....yn sono funzioni delle x allora la probabilita congiunta delley vale:

p(y1, y2, ....yn)dy1dy2...dyn = p(x1, x2, ...xn)|∂(x1, x2, ....xn)

∂(y1, y2, ....yn)|dy1, dy2, .....dyn ,

(3.33)dove |∂(......)/∂(....)| e il determinante Jacobiano delle x rispetto alley. Un esempio importante della formula suddetta e il metodo di Box-Muller per generare una distribuzione Gaussiana:

p(y)dy =1√2πe

−y22 dy . (3.34)

99

3.7. NUMERI RANDOM CAPITOLO 3. PROBABILITA E STATISTICA

Consideriamo adesso la trasformazione fra due distribuzioni uniformisu (0,1),x1,x2 e le due quantita y1,y2:

y1 =√−2ln x1cos(2πx2) , (3.35)

y2 =√−2ln x1sin(2πx2) . (3.36)

Oppure possiamo scrivere:

x1 = exp[−1

2(y21 + y22)] , (3.37)

x2 = exp1

2πarctan

y2y1

. (3.38)

Adesso il Jacobiano puo essere facilmente calcolato come:

∂(x1, x2)

∂(y1, y2)= −

[ 1

2πe−y

21/2]×[ 1

2πe−y

22/2]

. (3.39)

Abbiamo il prodotto delle funzioni y1 e y2 e quindi y e distribuita secon-do una distribuzione normale. Un ulteriore trucco consiste nel prenderev1 e v2 come ordinata ed ascissa di un punto random dentro il cerchiounitario intorno all’ origine. La somma dei loro quadrati R2 ≡ v21 +v22 e un generatore uniforme che puo essere usato per x1 mentre l’ an-golo che (v1, v2) definisce rispetto all’ asse v1 puo servire come angolorandom 2πx2. Il vantaggio e che seno e coseno possono essere scritticome v1/

√R2 e v2/

√R2 e quindi vengono eliminate le chiamate trigo-

nometriche. L’ implementazione pratica di questo algoritmo avvieneattraverso la FUNCTION GASDEV.

FUNCTION gasdev(idum)

INTEGER idum

REAL gasdev

CU USES ran1

INTEGER iset

REAL fac,gset,rsq,v1,v2,ran1

SAVE iset,gset

DATA iset/0/

if (iset.eq.0) then

1 v1=2.*ran1(idum)-1.

v2=2.*ran1(idum)-1.

100

CAPITOLO 3. PROBABILITA E STATISTICA 3.8. MOMENTI

rsq=v1**2+v2**2

if(rsq.ge.1..or.rsq.eq.0.)goto 1

fac=sqrt(-2.*log(rsq)/rsq)

gset=v1*fac

gasdev=v2*fac

iset=1

else

gasdev=gset

iset=0

endif

return

END

3.8 Momenti

Se gli elementi di un campione hanno la tendenza a raggrupparsi attornoad un valore particolare puo essere utile rappresentare l’ insieme dapochi numeri che sono i momenti; il piu importante e il valor medio,

x =1

N

N∑j=1

xj , (3.40)

che stima il valore centrale. Un altro importante valore e la varianzaempirica

s2 =1

N − 1

N∑j=1

(xj − x)2 , (3.41)

collegata ovviamente alla deviazione standard empirica

s =√s2 . (3.42)

In alcuni casi, tipo distribuzione di Lorentz, la deviazione standard nonesiste perche tende ad infinito; bisogna allora ricorrere ad uno stimatorepiu robusto che e la deviazione assoluta

ADEV (x1.......xN) =1

N

N∑j=1

|xj − x| . (3.43)

Questa definizione e stata poco adoperata dagli statistici per il fattoche essendo un valore assoluto i teoremi sono difficilmente provabili.

101

3.8. MOMENTI CAPITOLO 3. PROBABILITA E STATISTICA

Abbiamo poi i momenti di ordine superiore che sono in realta pocousati; sono in genere numeri adimensionali mentre ricordiamo che il va-lor medio e deviazione standard empirica sono dimensionali. Abbiamoquindi che la asimmetria e cosı definita

skew(x1....xN) =1

N

N∑j=1

[xj − xs

]3, (3.44)

dove s e la deviazione standard empirica. Un valore positivo dellaasimmetria significa che abbiamo una coda che si estende verso i va-lori positivi e negativo viceversa. Un altro parametro importante e lacurtosi

Kurt(x1....xN) ={ 1

N

N∑j=1

[xj − xs

]4}− 3 . (3.45)

Ricordiamo che questo parametro nel caso di una distribuzione Gaus-siana vale 0. Una maniera per calcolare la varianza che minimizza glierrori di arrotondamento e l’ algoritmo a due passi: prima il calcolo dix e poi quello di s2

s2(x1....xN) =1

N − 1

{ N∑j=1

(xj − x)2 − 1

N

[ N∑j=1

(xj − x)]2}

. (3.46)

La seconda somma e zero nel caso in cui x e esatto altrimenti corregge glierrori di arrotondamento del primo termine. I vari momenti del campio-ne che abbiamo teste definito sono implementati nella SUBROUTINEMOMENT

SUBROUTINE moment(data,n,ave,adev,sdev,var,skew,curt)

INTEGER n

REAL adev,ave,curt,sdev,skew,var,data(n)

INTEGER j

REAL p,s,ep

if(n.le.1)pause ’n must be at least 2 in moment’

s=0.

do 11 j=1,n

s=s+data(j)

11 continue

ave=s/n

adev=0.

102

CAPITOLO 3. PROBABILITA E STATISTICA 3.9. TEST DI STUDENT

var=0.

skew=0.

curt=0.

ep=0.

do 12 j=1,n

s=data(j)-ave

ep=ep+s

adev=adev+abs(s)

p=s*s

var=var+p

p=p*s

skew=skew+p

p=p*s

curt=curt+p

12 continue

adev=adev/n

var=(var-ep**2/n)/(n-1)

sdev=sqrt(var)

if(var.ne.0.)then

skew=skew/(n*sdev**3)

curt=curt/(n*var**2)-3.

else

pause ’no skew or kurtosis when zero variance in moment’

endif

return

END

3.9 Test di Student

Sovente ci troviamo di fronte due set di dati e vogliamo capire se ap-partengono o no alla stessa popolazione. Quando si pensa che due di-stribuzioni abbiano la stessa varianza empirica ma possibilmente mediediverse la variabile t di Student e calcolata in questa maniera

sD =

√∑iεA(xi − xA)2 +

∑iεB(xi − xB)2

NA +NB − 2

( 1

NA

+1

NB

), (3.47)

dove la sommatoria e estesa sui punti appartenenti ad un campione oda un altro e NA e NB sono il numero dei punti nel primo e secondo

103

3.9. TEST DI STUDENT CAPITOLO 3. PROBABILITA E STATISTICA

campione rispettivamente. Come secondo passaggio valutiamo t come:

t =xA − xBsD

. (3.48)

Terzo valutiamo la significativita di questo valore della t di Student conNA + NB - 2 gradi di liberta attraverso le due formule prima introdottee la SUBROUTINE BETAI (funzione beta incompleta ). La significa-tivita e un numero compreso fra 0 ed 1 e rappresenta la probabilita che|t| sia cosı grande per caso, nell’ ipotesi che i campioni abbiano la stessamedia. Quindi un valore piccolo della significativita ( 0.05 oppure 0.01) significa che la differenza osservata e molto significativa. Abbiamoquindi la SUBROUTINE TTEST dove data1 (1:n1 ) e data2 (1:n2 )sono i due vettori contenenti i dati dei due campioni.

SUBROUTINE ttest(data1,n1,data2,n2,t,prob)

INTEGER n1,n2

REAL prob,t,data1(n1),data2(n2)

CU USES avevar,betai

REAL ave1,ave2,df,var,var1,var2,betai

call avevar(data1,n1,ave1,var1)

call avevar(data2,n2,ave2,var2)

df=n1+n2-2

var=((n1-1)*var1+(n2-1)*var2)/df

t=(ave1-ave2)/sqrt(var*(1./n1+1./n2))

prob=betai(0.5*df,0.5,df/(df+t**2))

return

END

che fa uso della SUBROUTINE AVEVAR

SUBROUTINE avevar(data,n,ave,var)

INTEGER n

REAL ave,var,data(n)

INTEGER j

REAL s,ep

ave=0.0

do 11 j=1,n

ave=ave+data(j)

11 continue

ave=ave/n

104

CAPITOLO 3. PROBABILITA E STATISTICA 3.10. FIT CON RETTA

var=0.0

ep=0.0

do 12 j=1,n

s=data(j)-ave

ep=ep+s

var=var+s*s

12 continue

var=(var-ep**2/n)/(n-1)

return

END

3.10 Fit con retta

Supponiamo di avere una serie N di punti (xi,yi) e di volerli fittare conuna retta

y = a+ bx . (3.49)

La misura di quanto il modello aderisce alla realta e rappresentata dallafunzione chi-quadro:

χ2(a, b) =N∑j=1

(yi − a− bxisi

)2. (3.50)

Annullando le derivate prime rispetto ad a e b troviamo la condizio-ne di minimo del chi-quadro. Puo essere utile introdurre le sequentisommatorie:

S ≡N∑i=1

1

s2iSx ≡

N∑i=1

xis2i

Sy ≡N∑i=1

yis2i

(3.51)

Sxx ≡N∑i=1

x2is2i

Sxy ≡N∑i=1

xiyis2i

. (3.52)

Ricordando che il determinante dei coefficienti vale:

∆ ≡ SSxx − (Sx)2 , (3.53)

troviamo facilmente dalla teoria dei minimi quadrati che

a =Sy Sxx − Sx Sxy

∆b =

S Sxy − Sx Sy∆

. (3.54)

105

3.10. FIT CON RETTA CAPITOLO 3. PROBABILITA E STATISTICA

Ci restano ancora gli errori su a e b e dalla teoria della propagazionedegli errori otteniamo:

s2a = Sxx/∆ s2b = S/∆ . (3.55)

Una stima della bonta del fit e data dalla probabilita Q che il valoredel chi-quadro trovato sia casuale e

Q = gammq (N − 2

2,χ2

2) . (3.56)

Ricordiamo che GAMMQ e la nostra FUNCTION per la funzione gam-ma incompleta Q(a,x). Se q e piu grande di 0.1 il fit e accettabile. See piu grande di 0.001 il fit puo essere accettabile se gli errori non sonodistribuiti in maniera normale e sono stati sottostimati. Se q e minoredi 0.001 allora il modello o la procedura di stima degli errori sono sba-gliati. L’ implementazione pratica di queste formule avviene attraversola SUBROUTINE FIT dove in input abbiamo:ndata = numero di punti del campionex (1:ndata)= vettore delle xy (1:ndata)= vettore delle ysig (1:ndata)= vettore degli errori sulle ymwt = 0 allora applichiamo le formule senza errori

In output abbiamo invecea,b = parametri della rettasiga,sigb= errori sui parametri a e bchi2= chi-quadro finaleq= parametro di bonta del fit. Se mwt = 0 allora q=1

SUBROUTINE fit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)

INTEGER mwt,ndata

REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)

CU USES gammq

INTEGER i

REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,gammq

sx=0.

sy=0.

st2=0.

b=0.

if(mwt.ne.0) then

106

CAPITOLO 3. PROBABILITA E STATISTICA 3.10. FIT CON RETTA

ss=0.

do 11 i=1,ndata

wt=1./(sig(i)**2)

ss=ss+wt

sx=sx+x(i)*wt

sy=sy+y(i)*wt

11 continue

else

do 12 i=1,ndata

sx=sx+x(i)

sy=sy+y(i)

12 continue

ss=float(ndata)

endif

sxoss=sx/ss

if(mwt.ne.0) then

do 13 i=1,ndata

t=(x(i)-sxoss)/sig(i)

st2=st2+t*t

b=b+t*y(i)/sig(i)

13 continue

else

do 14 i=1,ndata

t=x(i)-sxoss

st2=st2+t*t

b=b+t*y(i)

14 continue

endif

b=b/st2

a=(sy-sx*b)/ss

siga=sqrt((1.+sx*sx/(ss*st2))/ss)

sigb=sqrt(1./st2)

chi2=0.

if(mwt.eq.0) then

do 15 i=1,ndata

chi2=chi2+(y(i)-a-b*x(i))**2

15 continue

q=1.

sigdat=sqrt(chi2/(ndata-2))

107

3.10. FIT CON RETTA CAPITOLO 3. PROBABILITA E STATISTICA

siga=siga*sigdat

sigb=sigb*sigdat

else

do 16 i=1,ndata

chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2

16 continue

q=gammq(0.5*(ndata-2),0.5*chi2)

endif

return

END

108

Capitolo 4

I frattali

In questi ultimi anni grazie all’ avvento dei calcolatori un capitolo chesembrava dimenticato della matematica ha avuto una evoluzione im-pressionante ed un numero molto alto di libri sui frattali ed affini estato scritto. Diciamo grosso modo che possiamo dividere in parti ilproblema:

1. Enunciazione del processo; spesso un algoritmo iterativo.

2. Sviluppo di uno pseudo-codice che indica le operazioni da svolgereindipendentemente dal linguaggio adoperato.

3. Traduzione pratica dello pseudo-codice in linguaggio di program-mazione

4. Uso di uno o piu pacchetti grafici che permettono l’ uscita finalesu terminale grafico, laser, etc...

Passiamo ora in rassegna alcuni di questi algoritmi seguendo unaspecie di ordine cronologico.

4.1 Henon mapping

All’ inizio degli anni 60 Michel Henon all’ Osservatorio di Nizza studian-do le orbite delle stelle in una galassia si rese conto che lievi modifichenei parametri iniziali portavano ad orbite irregolari. Prese quindi piedeil concetto di comportamento caotico nelle orbite stellari o planetarie.

109

4.1. HENON MAPPING CAPITOLO 4. I FRATTALI

Le curve di Henon non sono le classiche ellissi Newtoniane che gli astro-nomi usano da secoli. Nell’ universo vero l’ influenza di altri pianeti odi altre stelle portano a comportamenti poco prevedibili. Il modello diHenon che si sviluppa secondo il seguente mapping:

xn+1 = yn− ax2n (4.1)

yn+1 = bxn , (4.2)

con a = 1.4 e b = 0.3 , (4.3)

dice che le orbite planetarie possono essere non periodiche ed estre-mamente sensitive anche alla piu piccola influenza gravitazionale. Permaggiori dettagli sull’ attrattore di Henon consultare il libro [Ol 92] . Adifferenza dei modelli lineari che sembrano capaci di predirre il camminoesatto dei corpi celesti per l’ eternita, gli attrattori strani offrono la pos-sibilita di orbite stabili per molto tempo con un susseguente improvvisocambio di traiettoria vedi figura 4.1 . Riportiamo la SUBROUTINEHENON10 che implementa l’ attrattore di Henon.

SUBROUTINE HENON10 (xx,yy,max,kounter,n)

implicit REAL * 4 (A-H,O-Z)

dimension xx (1:max )

dimension yy (1:max )

common /quadro / xmin,ymin,delta

xmax=xmin + delta

ymax=ymin + delta

QA=1.4E+00

Qb=0.3E+00

J=0

k=0

18 continue

QZ=1.0-QA*(QX**2.0)+QY

QY=QB*QX

QX=QZ

k=k+1

IF ((QX.GT.xmin.AND.QX.LT.xmax).and.

& (QY.GT.ymin.AND.QY.LT.ymax))then

J=J+1

xx ( J ) = qx

yy ( J ) = qy

110

CAPITOLO 4. I FRATTALI 4.1. HENON MAPPING

Figura 4.1: Tipico plot di Henon

111

4.2. FEIGENBAUM MAPPING CAPITOLO 4. I FRATTALI

c fase di plot

end if

20 IF (K.LT.N) GO TO 18

kounter = j

return

END

4.2 Feigenbaum mapping

Questo mapping e collegato alla crescita della popolazione che assumia-mo essere x0 all’ inizio. La relazione ricorsiva gia adottata da Verhulste:

xn+1 = 4λxn(1− xn) , (4.4)

dove λ e un parametro di controllo che varia tra 0 ed 1. Questo map-ping al variare di λ incomincia con una soluzione, che poi diventanodue, che poi diventano quattro, etc... Questo comportamento prende ilnome di teoria delle biforcazioni, vedi figura 4.2. Il famoso matematicoFeigenbaum riesce poi anche a calcolare un numero che rappresenta latransizione al caos. Per maggiori dettagli sul mapping e sul numero diFeigenbaum consultare il libro [Ol 92] oppure il libro [Ba 88].

L’ implementazione pratica di questo algoritmo avviene attraversoil programma FEIGEN.

program feigen

c questo programma produce dei files di dati per calcolare lo splitting

c di Feigenbaum........................................................

implicit REAL (A-H,R-Z)

implicit real*8 (Q)

character * 1 mode,risposta

parameter ( nn = 200 )

parameter ( max= 1000 )

real* 8 fe ( 1 :max )

real* 4 fi ( 1 :max )

real* 4 fy ( 1:nn,1 : max )

integer *2 ify ( 1:nn )

real* 4 fx ( 1:nn )

mode =’M’

call init9 ( mode )

112

CAPITOLO 4. I FRATTALI 4.2. FEIGENBAUM MAPPING

Figura 4.2: Biforcazioni di Feigenbaum

113

4.2. FEIGENBAUM MAPPING CAPITOLO 4. I FRATTALI

go to 402

400 continue

type *,’Vuoi continuare Y/N ? ’

accept 1000,risposta

if ((risposta.eq.’N’).or.(risposta.eq.’n’)) then

go to 401

end if

1000 format (A1 )

402 continue

c qx0old ==> condizione iniziale ,ql0 ==> lambda iniziale

c qld ==> delta lamda ,npunti===>numero attrattori

c nn=divisioni in lamda

type *,’inserire limite inferiore + limite superiore ’

accept *,ql0,qxmax

qld = qxmax - ql0

qx0old=0.5

npunti=600

delta = 1.0 / max

qxmax=ql0+qld

qymin=0.0

qymax=1.0

do 111 ii=1,nn

qx0=qx0old

ql= ql0+ii*qld/nn

k=0

j=0

n=3000

maxfe=N-npunti

18 continue

QX=4*QL*QX0*(1.0-QX0)

qx0=qx

k=k+1

if (k.gt.maxfe) then

j=j+1

fe(j)=qx

end if

20 IF (K.LT.N) GO TO 18

k=0

do 67 j=1,npunti

114

CAPITOLO 4. I FRATTALI 4.3. MANDELBROT-SET

if (fe(j).lt.1.0D-10) then

fe(j)=1.0Q-10

end if

67 continue

call m01apf(fe,1,npunti,k)

fi(1)=fe(1)

mm=1

do 4 k=2,npunti

diff = abs (fe (k) - fi ( mm) )

if (diff .gt . delta) then

mm=mm+1

fi(mm)=sngl ( fe(k) )

endif

4 continue

do 3 j=1,mm

fy ( ii,j ) = fi ( j )

3 continue

ify( ii ) = mm

fx ( ii ) = sngl ( ql )

111 continue

xmin = sngl ( ql0 )

xmax = sngl ( ql0 + qld )

ymin = -0.1

ymax = 1.1

call feigengraf ( fy, fx ,nn ,max,xmin,xmax,ymin,ymax,ify)

go to 400

401 continue

call donepl

stop

END

4.3 Mandelbrot-set

Sia zn = xn + i yn un numero complesso definito dalla relazione ricorsiva

zn+1 = z2n + c , (4.5)

115

4.4. IMPLEMENTAZIONE DEI SET CAPITOLO 4. I FRATTALI

dove c =a + i b e un altro numero complesso. Se incominciamoil processo da z0 = 0 e e ripetiamo piu volte la relazione ricorsiva

notiamo che esistono due possibilita. La prima e che il valore assolutodi zn+1 esploda la crescere di n, la seconda e che rimanga finito a di-spetto di quante volte iteriamo la relazione. Il risultato finale dipendeovviamente dal valore di c usato. E quindi possibile dividere il set ditutti i numeri complessi {c} in due gruppi in accordo al loro andamentodopo un numero alto di iterazioni. L’ insieme di Mandelbrot e il gruppoche contiene tutti i c tali che

limx→∞| zn | = finito . (4.6)

4.4 Implementazione dei set

In genere il set di Mandelbrot e derivati vengono rappresentati tramiteil colore perche in tale maniera si originano delle immagini affascinantie ogni volta simili ma diverse nei particolari. Per maggiori dettaglisui vari set tipo Mandelbrot e Julia rimandiamo il lettore interessatoal libro di [PS88] . Per semplicita tratteremo solo la situazione sulmonitor di un computer. Supponiamo che il numero totale di pixelsdisponibile sia ( Nx×Ny ); dove Nx, ad esempio 640, sono i punti sullalinea orizzontale e Ny, ad esempio 480,sono i punti sulla linea verticaledi un tipico schermo di P.C. Possiamo invece anche lavorare su di unquadrato e allora il numero di pixels (quadratini) sara ( N ×N ) doveN e il numero considerato fissabile a piacimento; in questa maniera sieliminano eventuali errori di deformazione dell’ area. Identifichiamoadesso lo schermo con il piano complesso e ogni pixel o quadratino e unpunto del piano. Come al solito useremo la linea orizzontale per la partereale e la direzione verticale per la parte immaginaria. Per inserire unnumero complesso c =a + ib sullo schermo dobbiamo selezionare unorigine in basso a sinistra, un ∆ considerato ed il numero di pixels. Perfini pratici dobbiamo scegliere un valore nmax come il numero massimodi iterazioni che vogliamo eseguire; questo sara l’ infinito del nostroprocesso iterativo. Per generare il plot ogni numero complesso c esoggetto ad un processo ricorsivo partendo da z0 = 0. Alla fine di ogniiterazione il valore assoluto di zn e calcolato. Se | zn | > 2, il numerocomplesso c non fa parte del set di Mandelbrot, il processo iterativotermina e assegniamo un colore al pixel considerato secondo lo schemaseguente: sia Nesc il numero di iterazioni considerato ed Nc il numero

116

CAPITOLO 4. I FRATTALI 4.4. IMPLEMENTAZIONE DEI SET

di colori a disposizione, dove probabilmente Nc < Nmax; assoceremo alpixel il colore Ic secondo la seguente relazione

Ic = Nesc mod Nc , (4.7)

dove Ic e un numero nel range [0, Nc] ed ogni intero 0,1,2,.... corrispondead un particolare colore tipo nero, rosso,viola .... Se il valore di zn eancora minore di due il processo si ripete e se dopo Nmax iterazioni| zn | e ancora minore di due abbiamo un elemento appartenente al setdi Mandelbrot. Questo ultimo pixel puo essere colorato in nero: unatipica implementazione di questo set e riportata in figura 4.3 .

Riportiamo la subroutine UNO che implementa praticamente tra-mite tre do innestati il numero Nesc ( nel nostro caso NDIM ) per ognipixel del piano.

subroutine uno ( matrix,traspo,itraspo)

implicit real* 4 (a-h,o-z)

integer * 2 ndim,i,j,iterations

integer * 4 im

complex * 8 z,imm,complex

parameter ( pixels = 600 )

parameter ( nvector = pixels * pixels)

dimension matrix (1:nvector)

dimension traspo (1:10)

dimension itraspo (1:10)

parameter(imm=(0,1))

c coordinate di partenza e valore del lato

xmin = traspo (1)

ymin = traspo (2)

delta = traspo (3)

npix = itraspo(1)

iterations = itraspo (2)

npix1 = npix-1

xmax = xmin + delta

ymax = ymin + delta

c griglia di lavoro + numero di pixels

im = 0

do i= 0,npix1

cy=ymin + i*delta/npix

do j=0,npix1

117

4.4. IMPLEMENTAZIONE DEI SET CAPITOLO 4. I FRATTALI

Figura 4.3: Tipico set di Mandelbrot

118

CAPITOLO 4. I FRATTALI 4.4. IMPLEMENTAZIONE DEI SET

cx=xmin+j*delta/npix

complex = cx + imm*cy

c dimensione numero complesso

z=0

ndim = 0

do while (cabs(z) .le. 2.0 .AND. ndim .LT. iterations)

z=z*z+complex

ndim=ndim+1

end do

im = im +1

matrix (im) = ndim

end do

end do

return

end

Il colore viene invece calcolato dalla subroutine DUETUTTO.

subroutine duetutto( matrix,matrix8,traspo,itraspo)

implicit real * 4 (a-h,o-z)

parameter ( npixels=800 )

parameter ( nvector = npixels*npixels)

dimension matrix (1:nvector )

dimension matrix8 (1:nvector )

dimension traspo (1:10)

dimension itraspo (1:10)

a = traspo (1)

npix = itraspo(1)

ncol = itraspo(3)

npix2 = npix *npix

do j = 1 , npix2

ndim = matrix ( j)

ndim = ndim -1

if (ndim.lt.64) then

matrix8 (j) = ndim

elseif (ndim.ge.64 ) then

icolor = mod ( ndim,64 )

matrix8 (j) = icolor

end if

end do

119

4.5. SET DI JULIA E DERIVATI CAPITOLO 4. I FRATTALI

return

end

4.5 Set di Julia e derivati

Gli insiemi di Julia si ottengono dalla stessa relazione ricorsiva del set diMandelbrot e si tratta di selezionare dei punti alla frontiera del mede-simo; questi punti dell’ altro set danno origine a successioni periodichedalle forme strane e bizzarre. I pixels sullo schermo del computer rap-presentano ora i punti nel piano complesso per tutti i possibili valoridi z0 ed ogni insieme e caratterizzato da un valore di c. Riportia-mo nella figura 4.4 il set di Julia corrispondente a c = -0.74543 + i0.11301. Un’ altra relazione ricorsiva cugina di quella di Mandelbrot,che possiamo chiamare M2, e

zn+1 = z3n + c , (4.8)

che possiamo chiamare M3. Essa origina strutture leggermente diversedal suo parente di grado minore e riportiamo nella figura 4.5 un tipicoplot. Possiamo poi passare in maniera continua da M2 a M3 attraversola seguente relazione ricorsiva

zn = tz3 + (1− t)z2 + c . (4.9)

Dove il parametro t varia con continuita fra 0 ed 1.Sembra che alcuni valori di t siano critici e potete partire con t

=0.12595. Riportiamo in figura 4.6 un tipico plot di di questo setregolabile.

4.6 Le trasformazioni affini

Il metodo IFS ( Iterated Function System ) permette di produrre im-magini di curve frattali, foglie, nuvole, etc...

Ricordiamo che una trasformazione affine wi e data da

wi = [xy ] = [ a bc d ] [xy ]+[ ef ] = [ ax+ by+ ecx+ dy+ f ] .(4.10)

120

CAPITOLO 4. I FRATTALI 4.6. LE TRASFORMAZIONI AFFINI

Figura 4.4: Tipico plot di Julia

121

4.6. LE TRASFORMAZIONI AFFINI CAPITOLO 4. I FRATTALI

Figura 4.5: Tipico plot M**3

122

CAPITOLO 4. I FRATTALI 4.6. LE TRASFORMAZIONI AFFINI

Figura 4.6: Tipico plot regolabile

123

4.6. LE TRASFORMAZIONI AFFINI CAPITOLO 4. I FRATTALI

Tabella 4.1: Tabella dei coefficienti

w a b c d e f p

1 0 0 0 0.16 0 0 0.012 0.85 0.04 −0.04 −0.85 0 1.6 0.853 0.20 −0.26 0.23 0.22 0 1.6 0.074 −0.15 0.28 0.26 0.24 0 0.44 0.07

Inoltre ogni trasformazione e caratterizzata da una certa probabilitadi avvenire pi associata a wi. Ovviamente essendo probabilita dovremoavere:

p1 + p2 + .....pn = 1 e pi > 0 per i = 1, 2, ...n . (4.11)

Riportiamo in tabella i valori dei coefficienti che permettono di ottenereuna figura caratteristica.

Per capire le componenti di queste strutture complicate possiamostudiare il significato geometrico dei coefficienti a,b,c,d,e ed f. Per par-tire osserviamo che e ed f sono semplici traslazioni nelle direzioni x edy. Adesso scriviamo a = r × cos ( θ ), b = -s × sen ( ψ ), c = r×sen ( θ) e d = s × cos (ψ ). E facile realizzare che θ e l’ angolo delquale ruotiamo l’ asse x e ψ e l’ angolo del quale ruotiamo y, r e fattoreper il quale sono moltiplicate le distanze lungo x ed s e il fattore per ilquale moltiplicare le distanze lungo y. Prendendo θ = ψ abbiamo unariflessione attravesrso l’ asse x e prendendo s negativo si origina unariflessione attraverso l’ asse y. Possiamo quindi riesprimere i parame-tri a, b, c,d attraverso r, s, θ e ψ; continuando possiamo mettere h ek al posto di e ed f. Per maggiori dettagli sulle trasformazioni affiniconsultare il libro di [Ba 88] . Riportiamo nella figura 4.7 seguente unatipica implementazione pratica. L’ implementazione pratica avvieneattraverso il programma FELCI.

program felci

implicit real *8 (a-h,o-z)

logical * 1 prob1,prob2,prob3,prob4

real * 4 xx , yy

parameter ( nvec = 150 000 )

124

CAPITOLO 4. I FRATTALI 4.6. LE TRASFORMAZIONI AFFINI

Figura 4.7: La felce

125

4.6. LE TRASFORMAZIONI AFFINI CAPITOLO 4. I FRATTALI

parameter ( minimum = 100 )

dimension xx (1 : nvec )

dimension yy (1: nvec )

integer g05dyf

call g05cbf ( 0 )

type *,’numero di punti ?’

accept *,npoints

x = 0.0

y = 0.0

do j = 1, npoints

intero = g05dyf ( 1, 100 )

prob1 = intero.eq.1

prob2 = (intero.gt.1 ).and.(intero.le.86 )

prob3 = (intero.gt.86).and.(intero.le.93 )

prob4 = (intero.gt.93).and.(intero.le.100)

if ( prob1 ) then

itipo = 1

call affine ( itipo,x,y,xnew,ynew )

elseif ( prob2 ) then

itipo = 2

call affine ( itipo,x,y,xnew,ynew )

elseif ( prob3 ) then

itipo = 3

call affine ( itipo,x,y,xnew,ynew )

elseif ( prob4 ) then

itipo = 4

call affine ( itipo,x,y,xnew,ynew )

end if

if ( j.gt. minimum ) then

jvec = jvec + 1

xx ( jvec ) = sngl ( xnew )

yy ( jvec ) = sngl ( ynew )

end if

x = xnew

y = ynew

end do

call affinegraf ( xx, yy , jvec ,npoints,minimum,nvec)

stop

end

126

CAPITOLO 4. I FRATTALI 4.6. LE TRASFORMAZIONI AFFINI

subroutine affine ( itipo,x,y,xnew,ynew )

implicit real *8 ( a-h,o-z)

if (itipo.eq.1) then

a = 0.

b = 0.

c = 0.

d = 0.16

e = 0.

f = 0.

elseif (itipo.eq.2) then

a = 0.85

b = 0.04

c = -0.04

d = 0.85

e = 0.0

f = 1.6

elseif (itipo.eq.3) then

a = 0.2

b = -0.26

c = 0.23

d = 0.22

e = 0.0

f = 1.6

elseif (itipo.eq.4) then

a = -0.15

b = 0.28

c = 0.26

d = 0.24

e = 0.00

f = 0.44

end if

xnew = a * x + b * y + e

ynew = c * x + d * y + f

return

end

127

4.7. IL SISTEMA L CAPITOLO 4. I FRATTALI

4.7 Il sistema L

Questo sistema e stato inventato dal botanico tedesco Lindenmayer (vedi il libro di [PL 90]) per spiegare teoricamente le strutture ricorrentidelle piante. Ai fini pratici possiamo sintetizzare l’ algoritmo medianteil moto di una tartaruga. Lo stato della tartaruga e definito da unatripletta ( x,y, α ), dove le coordinate cartesiane ( x, y ) rappresentanola posizione della tartaruga e l’ angolo α indica la direzione di motodella tartaruga. Dato l’ incremento d di ogni step ( possiamo porlouguale ad uno ) e l’ incremento di angolo δ la tartaruga risponde aiseguenti comandi

• F Si muove avanti di uno step d. Lo stato della tartaruga diventa( x’, y’, α ) dove x’ = x + d cosα e y’ = y + d sin α. Una lineaviene tracciata fra i punti (x,y) e (x’,y’).

• f Si muove avanti di una lunghezza d senza tracciare una linea.

• + Si gira a sinistra di un angolo δ. Il prossimo stato della tartaru-ga e (x,y, α + δ ) e quindi gli angoli crescono in senso antiorario.

• - Ruota a sinistra di un angolo δ. Il prossimo stato della tartarugae (x,y,α - δ ).

• [ Congela lo stato della tartaruga. Le informazioni correnti, po-sizione e angolo sono memorizzate per poi essere riesumate piuavanti.

• ] Scongela lo stato della tartaruga. Le informazioni vengono riesu-mate e diventano quelle attuali della tartaruga. Non viene tiratanessuna linea e quasi sempre la tartaruga cambia posizione.

Data una stringa iniziale w, lo stato iniziale della tartaruga con (x0, , y0,α0 ) e fissati i parametri d e δ, tramite una o piu regole diproduzione otteniamo la posizione successiva e cosı via. Il problema equindi di ideare una stringa iniziale e regole di trasformazione tali daottenere le curve volute. Riportiamo come esempio una pianta curvatadal vento, vedi figura 4.8 .n = 4δ = 22.5w = F

128

CAPITOLO 4. I FRATTALI 4.7. IL SISTEMA L

Figura 4.8: Pianta e vento di tramontana

129

4.7. IL SISTEMA L CAPITOLO 4. I FRATTALI

F → FF-[-F+F+F]+[+F-F-F]

A questo punto ricordiamo che sono attualmente disponibili unaquarantina di curve fra frattali e piante. Questa particolare stringaviene implementata attraverso la SUBROUTINE MAKEPIANTA3.

subroutine makepianta3 (stringa,max,it,ifine)

implicit real * 4 (a-h,o-z)

parameter (new = 200 000 )

character * 1 stringa ,stringanew

dimension stringa (1:max)

dimension stringanew (1:new)

ifine = 1

stringa (1) = ’F’

do jit =1,it

nuovo = 0

do j = 1 ,ifine

if (stringa (j).eq.’F’) then

nuovo = nuovo +1

stringanew (nuovo) =’F’

nuovo= nuovo +1

stringanew (nuovo) =’F’

nuovo= nuovo +1

stringanew (nuovo) =’+’

nuovo= nuovo +1

stringanew (nuovo) =’[’

nuovo= nuovo +1

stringanew (nuovo) =’+’

nuovo= nuovo +1

stringanew (nuovo) =’F’

nuovo = nuovo +1

stringanew (nuovo )=’-’

nuovo= nuovo +1

stringanew (nuovo )=’F’

nuovo = nuovo +1

stringanew (nuovo )=’-’

nuovo= nuovo +1

stringanew (nuovo )=’F’

nuovo= nuovo +1

stringanew (nuovo )=’]’

130

CAPITOLO 4. I FRATTALI 4.7. IL SISTEMA L

nuovo= nuovo +1

stringanew (nuovo )=’-’

nuovo = nuovo +1

stringanew (nuovo) =’[’

nuovo= nuovo +1

stringanew (nuovo) =’-’

nuovo= nuovo +1

stringanew (nuovo) =’F’

nuovo= nuovo +1

stringanew (nuovo) =’+’

nuovo= nuovo +1

stringanew (nuovo) =’F’

nuovo= nuovo +1

stringanew (nuovo) =’+’

nuovo= nuovo +1

stringanew (nuovo )=’F’

nuovo= nuovo +1

stringanew (nuovo )=’]’

else

nuovo = nuovo +1

stringanew (nuovo) = stringa(j)

end if

end do

do j= 1,nuovo

stringa(j) = stringanew (j)

end do

ifine = nuovo

end do

return

end

Invece la stringa finale risultante viene decodificata attraverso la SUB-ROUTINE GX.

subroutine gx(stringa,max,ifine,agl,xmax,xmin,ymax,ymin,ae)

implicit real * 4 (a-h,o-z)

character * 1 stringa,car

dimension stringa (1:max)

dimension stackx ( 1:1000)

dimension stacky ( 1:1000)

131

4.7. IL SISTEMA L CAPITOLO 4. I FRATTALI

dimension stackae ( 1:1000)

xmax = - 10.**20

xmin = 10.**20

ymax = - 10.**20

ymin = 10.**20

x = 0.

y = 0.

numero = 0

do j =1,ifine

if (x.gt.xmax ) then

xmax = x

end if

if (x.lt.xmin ) then

xmin = x

end if

if (y.gt.ymax ) then

ymax = y

end if

if (y.lt.ymin ) then

ymin = y

end if

car = stringa (j)

cccccccccccc inizio sostituzioni

if (car.eq.’F’) then

xnew = x + cosd ( ae )

ynew = y + sind ( ae )

x = xnew

y = ynew

end if

if (car.eq.’f’) then

xnew = x + cosd ( ae )

ynew = y + sind ( ae )

x = xnew

y = ynew

end if

if (x.gt.xmax ) then

xmax = x

end if

if (x.lt.xmin ) then

132

CAPITOLO 4. I FRATTALI 4.7. IL SISTEMA L

xmin = x

end if

if (y.gt.ymax ) then

ymax = y

end if

if (y.lt.ymin ) then

ymin = y

end if

if (car.eq.’+’) then

ae = ae - agl

end if

if (car.eq.’-’) then

ae = ae + agl

end if

if (car.eq.’[’) then

numero = numero + 1

stackx ( numero ) = x

stacky ( numero ) = y

stackae (numero ) = ae

end if

if (car.eq.’]’) then

x = stackx ( numero )

y = stacky ( numero )

ae = stackae ( numero )

numero = numero - 1

end if

ccccccccccccccc fine sostituzioni ccccccccccccccccccccccccccc

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

end do

return

end

133

4.7. IL SISTEMA L CAPITOLO 4. I FRATTALI

134

Capitolo 5

Sistemi Monte Carlo

In questi ultimi anni abbiamo assistito ad un crescente numero di si-mulazioni di sistemi fisici mediante modelli al calcolatore. Il campo eenorme ed in continua evoluzione e ci limiteremo ad una scelta ran-dom di alcune topiche; ah non dimentichiamoci che il nome deriva dalfamoso casino...

5.1 DLA

Il nome prende le iniziali da Diffusion Limited Aggregation ed e un ti-pico esempio di crescita/aggregazione che simula il comportamento dicolonie di batteri, esperimenti fra fluidi con viscosita diverse e scaricheelettriche nei gas. La superficie considerata viene identificata con unlattice consistente in N×N pixels. All’ inizio il nostro aggregato consi-ste in una sola particella al centro della griglia. La diffusione delle altreparticelle viene simulata da una particella che entra da uno dei quattrolati a seconda del numero random R1[1, 4]. La particella si muove poidi moto random sul lattice: questo significa che ad ogni step ha quattropossibili movimenti scelti con un processo random R2[1, 4]. Le quattrodirezioni di movimento sono:

1. destra,

2. alto,

3. sinistra,

4. basso.

135

5.1. DLA CAPITOLO 5. SISTEMI MONTE CARLO

Ogni volta che viene raggiunta una locazione nuova viene controllatose uno dei quattro vicini ( di nuovo destra, alto, sinistra e basso ) eoccupato. Se risulta occupato l’ aggregato aumenta di una unita e ilcammino random ricomincia con un altra particella su di un altro latodel quadrato. Se invece il punto occupato non ha vicini appartenentiall’ aggregato il cammino random continua; il processo si ferma quandoviene raggiunto un punto sulla frontiera e ricomincia ovviamente con unaltra particella su di un altro lato. Si possono poi introdurre ulterioriregole che permettono di diminuire il tempo di calcolo; in particolarepossiamo ingrandire progressivamente il lato del quadrato interno dalquale partiamo e idem dicasi per un cuscinetto detto killing zone piugrande del quadrato delle partenze superato il quale il processo cominciadi nuovo senza aspettare di raggiungere le frontiere naturali del lattice.Per maggiori dettagli sulla teoria della DLA consultare il libro [Vi 92].Riportiamo nella figura 5.1 un tipico esempio di crescita su di un latticecon 500 pixels.

Questo algoritmo viene implementato attraverso la SUBROUTINEWO:

subroutine wo(griglia ,ndim ,l_dist,k_dist,m_dist,n_max,n_vero)

implicit real * 4 ( a-h,o-z)

logical * 1 griglia , occupato

external g05dyfzan

integer g05dyfzan

dimension griglia(ndim,ndim)

dimension ivicinox ( 1 : 4 )

dimension ivicinoy ( 1 : 4 )

dimension iarrivo ( 1 : 4 )

c zero the lattice and mark the center as occupied

do i=1,ndim

do j=1,ndim

griglia(j,i)=.false.

end do

end do

npts=1

mid=ndim/2

griglia(mid,mid)=.true.

c define initial launching and kill boundaries

lx_min=mid-l_dist

136

CAPITOLO 5. SISTEMI MONTE CARLO 5.1. DLA

Figura 5.1: Tipico esempio di DLA

137

5.1. DLA CAPITOLO 5. SISTEMI MONTE CARLO

lx_max=mid+l_dist

ly_min=mid-l_dist

ly_max=mid+l_dist

kx_min=lx_min-k_dist

kx_max=lx_max+k_dist

ky_min=ly_min-k_dist

ky_max=ly_max+k_dist

c set an artificial occupation boundary at the start

mx_min=mid-m_dist

mx_max=mid+m_dist

my_min=mid-m_dist

my_max=mid+m_dist

cc

cccc abbiamo fissato le zone di lancio

cc e quelle di killing

idum = 33333

kgen = 77777

kdir = 55555

nn = 1

ilato = 0

itotale= 0

ccc inizio zona di lancio

40 continue

ilato = g05dyfzan ( 1 , 4 ,idum )

itotale = itotale + 1

c ilato = ilato + 1

c if (ilato.gt.4) then

c ilato = 1

c end if

c

c lx_min,ly_max 2 lx_max,ly_max

c -------------------------

c | |

c | |

c | zona di |

c 1 | | 3

c | lancio |

c | |

138

CAPITOLO 5. SISTEMI MONTE CARLO 5.1. DLA

c -------------------------

c lx_min,ly_min 4 lx_max,ly_min

c

if ( ilato.eq. 1 ) then

ix = lx_min

iy = g05dyfzan ( ly_min,ly_max,kgen )

elseif ( ilato.eq. 2 ) then

ix = g05dyfzan ( lx_min,lx_max,kgen )

iy = ly_max

elseif ( ilato.eq. 3 ) then

ix = lx_max

iy = g05dyfzan ( ly_min,ly_max,kgen )

elseif ( ilato.eq. 4 ) then

ix = g05dyfzan ( lx_min,lx_max,kgen )

iy = ly_min

end if

c adesso impostiamo la direzione di moto

100 continue

idir = g05dyfzan ( 1 , 4 , kdir )

c

c

c 3

c |

c |

c |

c 4 ------ * ------ 2

c |

c |

c |

c 1

c

if (idir .eq . 1 ) then ! si muove giu‘ in basso

ix = ix

iy = iy - 1

elseif ( idir.eq.2 ) then ! si muove a destra

ix = ix +1

iy = iy

139

5.1. DLA CAPITOLO 5. SISTEMI MONTE CARLO

elseif ( idir.eq.3 ) then ! si muove in alto

ix = ix

iy = iy + 1

elseif ( idir.eq.4 ) then ! si muove a sinistra

ix = ix -1

iy = iy

end if

c

cc adesso ci siamo mossi !!!!!

c

c

c adesso facciamo un controllo sulla killing zone

c se abbiamo raggiunto la frontiera ricominciamo

c dalla zona di lancio

if ((iy.eq.ky_min).or.(iy.eq.ky_max)) then ! limite alto

! e basso

go to 40

elseif ((ix.eq.kx_min).or.(ix.eq.kx_max))then ! limite destro

! e sinistro

go to 40

end if

c

cc costruzione dei quattro vicini !

c

c -------------

c

c

c

c

c

c

ivicinox (1) = ix

ivicinoy (1) = iy - 1

c primo vicino

ivicinox (2) = ix + 1

140

CAPITOLO 5. SISTEMI MONTE CARLO 5.1. DLA

ivicinoy (2) = iy

c secondo vicino

ivicinox (3) = ix

ivicinoy (3) = iy +1

c terzo vicino

ivicinox (4) = ix - 1

ivicinoy (4) = iy

c quarto vicino

c

cc adesso vediamo se uno dei quattro vicini

c e‘ occupato

occupato =.false.

do jj = 1, 4

if ( griglia ( ivicinox ( jj), ivicinoy (jj)) ) then

if ( .not.griglia (ix,iy)) then

nn = nn + 1

if ( nn.eq.( n_max+1 ) ) then

nn = nn -1

type *,’Abbiamo raggiunto n_max’

go to 950

end if

griglia (ix,iy) = . true .

occupato = . true .

iarrivo (ilato )= iarrivo ( ilato ) +1

c

ccc ridefinizione dei massimi e dei minimi in x

c

if (ix.gt.mx_max) then

mx_max = ix

end if

if (ix.lt.mx_min) then

mx_min = ix

end if

c

ccc ridefinizione dei massimi e dei minimi in y

c

if (iy.gt.my_max) then

141

5.1. DLA CAPITOLO 5. SISTEMI MONTE CARLO

my_max = iy

end if

if (iy.lt.my_min) then

my_min = iy

end if

end if

end if

c

ccc ridefinizione eventuale della zona di lancio

c

if (( mx_min -lx_min) .eq. 1 ) then

go to 200

elseif (( lx_max - mx_max ) .eq. 1 ) then

go to 200

elseif (( ly_max - my_max ) .eq. 1 ) then

go to 200

elseif (( my_min -ly_min ) .eq. 1 ) then

go to 200

end if

go to 300

200 continue

c

ccccccc

c

lx_min = lx_min - l_dist

lx_max = lx_max + l_dist

ly_min = ly_min - l_dist

ly_max = ly_max + l_dist

idum = idum +1

kgen = kgen +1

ddir = kdir + 1

c

cccc ma potremmo avere raggiunto la frontiera

c

c

ccccccc

c

kx_min=lx_min-k_dist

142

CAPITOLO 5. SISTEMI MONTE CARLO 5.1. DLA

kx_max=lx_max+k_dist

ky_min=ly_min-k_dist

ky_max=ly_max+k_dist

cc

ccccccccc ma potremmo avere raggiunto le frontiere

cc

if (kx_min.le. 1 ) then

type *,’Abbiamo raggiunto la frontiera’

go to 950

elseif (ky_min.le. 1 ) then

type *,’Abbiamo raggiunto la frontiera’

go to 950

elseif (ky_max.ge. ndim ) then

type *,’Abbiamo raggiunto la frontiera’

go to 950

elseif (kx_max.ge. ndim ) then

type *,’Abbiamo raggiunto la frontiera’

go to 950

end if

300 continue

end do

if (.not.occupato ) then

go to 100

end if

go to 40

950 continue

call esci ( nn,n_vero)

isum = 0

do j = 1, 4

type *,’numero di arrivi dal lato ’,j,’ =’,iarrivo (j)

isum = isum + iarrivo (j)

end do

fraz = isum * 1.0 / itotale

type *,’frazione arrivata sui totali partiti =’,fraz

return

end

143

5.2. AGGREGAZIONE BALLISTICACAPITOLO 5. SISTEMI MONTE CARLO

5.2 Aggregazione ballistica

Questo processo di aggregazione cerca di spiegare fenomeni in cui entrain gioco la gravita, tipo deposizione di zinco su di un catodo, esperi-menti di soluzione chimica, depositi di ossido di manganese su calcari,neve che cade su di un vetro. Come al solito lavoriamo su di un latticecomposto da N ×N pixels e le particelle partono da un punto randomR1 [1,N] sul lato in alto. Inizialmente l’ aggregato consiste in una par-ticella posta sul lato in basso in posizione centrale ( deposito a cono )oppure in piu di una particella ( deposito ad alberi ). Ad ogni sposta-mento delle particelle in caduta vengono controllati i quattro vicini. Seun vicino e occupato l’ aggregato aumenta di una unita e il processoriparte con un’ altra particella sul lato in alto. Se invece nessuno deiquattro vicini e occupato il processo riparte con il moto random. Sevengono raggiunti i lati-frontiera sinistro oppure destro oppure inferio-re il processo si ferma e riparte dall’ alto. Anche quı per velocizzareil tempo di esecuzione del programma le particelle partono sul lato su-periore variabile. Riportiamo nelle figure 5.2 e 5.3 due tipici esempi diaggregazione ballistica.

L’ implementazione pratica avviene attraverso la SUBROUTINEBA.

subroutine ba (occ ,ndim ,l_dist,k_dist,m_dist,n_max,n_vero,nba)

implicit real * 4 ( a-h,o-z)

logical * 1 occ , occupato

external g05dyfzan

integer g05dyfzan

dimension occ(ndim,ndim)

dimension ivicinox ( 1 : 4 )

dimension ivicinoy ( 1 : 4 )

c zero the lattice and mark the center as occupied

do i=1,ndim

do j=1,ndim

occ(j,i)=.false.

end do

end do

idum = 33333

do i =1 , nba

ipunto = g05dyfzan ( 1 , ndim ,idum )

144

CAPITOLO 5. SISTEMI MONTE CARLO5.2. AGGREGAZIONE BALLISTICA

Figura 5.2: Deposito ballistico con un punto iniziale

145

5.2. AGGREGAZIONE BALLISTICACAPITOLO 5. SISTEMI MONTE CARLO

Figura 5.3: Deposito ballistico con piu punti iniziali

146

CAPITOLO 5. SISTEMI MONTE CARLO5.2. AGGREGAZIONE BALLISTICA

if (nba.eq.1) then

ipunto = ndim / 2

end if

occ (ipunto,1) =.true.

end do

c define initial launching and kill boundaries

lx_min = 1

lx_max = ndim

ly_max= 1 + l_dist

ly_min= 1

ky_min = 1

ky_max= ly_max+k_dist

kx_max = ndim

kx_min = 1

c set an artificial occupation boundary at the start

my_max= 1 + m_dist

cc

cccc abbiamo fissato le zone di lancio

cc e quelle di killing

idum = 33333

kgen = 77777

kdir = 55555

nn = 1

ccc inizio zona di lancio

40 continue

c

ccccccccc

c lx_min,ly_max ..... lx_max,ly_max

c -------------------------

c | |

c | |

c | zona di |

c | |

c | lancio |

c | |

c

ix = g05dyfzan ( lx_min,lx_max,kgen )

147

5.2. AGGREGAZIONE BALLISTICACAPITOLO 5. SISTEMI MONTE CARLO

iy = ly_max

c adesso impostiamo la direzione di moto

100 continue

idir = g05dyfzan ( 1 , 3 , kdir )

c

c 3 ------ * ------ 2

c |

c |

c |

c 1

c

if (idir .eq . 1 ) then ! si muove giu‘ in basso

ix = ix

iy = iy - 1

elseif ( idir.eq.2 ) then ! si muove a destra

ix = ix +1

iy = iy

elseif ( idir.eq.3 ) then ! si muove a sinistra

ix = ix -1

iy = iy

end if

c

cc adesso ci siamo mossi !!!!!

c

c

c adesso facciamo un controllo sulla killing zone

c se abbiamo raggiunto la frontiera ricominciamo

c dalla zona di lancio

if (iy.eq.ky_min) then ! limite alto

! e basso

go to 40

elseif ((ix.eq.kx_min).or.(ix.eq.kx_max))then ! limite destro

! e sinistro

go to 40

148

CAPITOLO 5. SISTEMI MONTE CARLO5.2. AGGREGAZIONE BALLISTICA

end if

c

cc costruzione dei tre vicini !

c

c -------------

c

c

c

c

c

c

ivicinox (1) = ix

ivicinoy (1) = iy - 1

c primo vicino

ivicinox (2) = ix + 1

ivicinoy (2) = iy

c secondo vicino

ivicinox (3) = ix - 1

ivicinoy (3) = iy

c terzo vicino

c

cc adesso vediamo se uno dei quattro vicini

c e‘ occupato

occupato =.false.

do jj = 1, 3

if ( occ ( ivicinox ( jj), ivicinoy (jj)) ) then

if ( .not.occ (ix,iy)) then

nn = nn + 1

if ( nn.eq.( n_max+1 ) ) then

nn = nn -1

call esci ( nn,n_vero)

type *,’Abbiamo raggiunto n_max’

return

end if

occ (ix,iy) = . true .

occupato = . true .

c

149

5.3. PERCOLAZIONE CAPITOLO 5. SISTEMI MONTE CARLO

ccc ridefinizione del massimo in y

c

if (iy.gt.my_max) then

my_max = iy

ly_max = my_max + l_dist

end if

if (ly_max.ge. ndim ) then

call esci ( nn,n_vero)

type *,’Abbiamo raggiunto la frontiera’

return

end if

end if

end if

end do

if (.not.occupato ) then

go to 100

end if

go to 40

end

5.3 Percolazione

Anche quı grazie all’ avvento dei computers e stato aperto un nuovocapitolo al confine fra fisica e matematica. Prendiamo in esame unlattice quadrato con N ×N occupazioni disponibili e introduciamo unaprobabilita di occupazione dei punti del lattice p varianti fra 0 ed 1.Questo significa che ad ogni p avremo un numero M di punti occupatisul lattice pari a:

M = NINT (p×N ×N) . (5.1)

Possiamo generare queste M occupazioni con due processi random{Ix = R1[1, N ]Iy = R2[1, N ]

, (5.2)

controllando ogni volta che il punto del lattice ( Ix, Iy) non sia giaoccupato. Alla fine del processo avremo M punti occupati sul lattice (N×N ). Notiamo che quando P ≥ 0.59 e possibile trovare un cammino

150

CAPITOLO 5. SISTEMI MONTE CARLO 5.3. PERCOLAZIONE

di siti occupati che connette alto e basso e destra e sinistra. Diciamo aquesto punto che il cluster percola attraverso il sistema come l’ acquapercola attraverso una caffettiera. Questa probabilita prende il nomedi probabilita critica pc e gioca un ruolo cruciale in tutta la teoria.Riportiamo nella figura 5.4 un tipico esempio di occupazione randomdi un cluster di 300 pixels con p= 0.1 .

Le applicazioni pratiche della percolazione possono essere gli in-cendi nelle foreste, le transizioni di fase tipo ferro-magnetismo ↔ pa-ramagnetismo, conduttivita, etc. ; per maggiori dettagli consultateil libro [TF 92]. Praticamente possiamo implementare l’ algoritmomediante la SUBROUTINE FAI-PUNTI.

subroutine fai_punti (kubo,mmm,perc)

implicit real * 4 ( a-h ,o-z)

logical * 1 kubo

dimension kubo ( 0 : mmm , 0 : mmm )

external g05dyfzan

integer g05dyfzan

do ix = 0, mmm

do iy = 0, mmm

kubo ( ix, iy ) = .false.

end do

end do

itotal = ( mmm + 1 ) *( mmm + 1)

iselected = nint ( perc * itotal )

icambiati = 0

idum = 777777

do jj = 1 , iselected

50 continue

ix = g05dyfzan (0,mmm,idum)

iy = g05dyfzan (0,mmm,idum)

if ( kubo ( ix,iy )) then

icambiati = icambiati + 1

go to 50

end if

kubo (ix,iy) =.true.

end do

return

end

151

5.3. PERCOLAZIONE CAPITOLO 5. SISTEMI MONTE CARLO

Figura 5.4: Percolazione con p=0.1

152

CAPITOLO 5. SISTEMI MONTE CARLO 5.3. PERCOLAZIONE

Il problema di trovare il cluster di connessione puo essere cosı risolto.Prendiamo un punto occupato del lattice, eliminiamolo dal lattice emettiamolo in un vettore. Abbiamo adesso quattro vicini e controlliamoquali sono quelli occupati: essi finiranno nel vettore. Calcoliamo adessoi vicini occupati da tutti gli elementi del vettore e cosı via. Questo eun cluster e il processo si ferma quando non troviamo piu nessun vicinooccupato. Il cluster di connessione viene trovato quando quattro puntidi esso appartengono ai lati sinistro, destro alto e basso; vedi figura 5.5.

Questo algoritmo di ricerca del cluster viene implementato attra-verso la SUBROUTINE FAICLUSTER.

subroutine fai_cluster ( kubo, mmm ,si ,cluster)

implicit real * 4 (a-h,o-z)

integer *4 ifine

logical *1 kubo,condizionea,condizioneb

logical *1 condizionex

logical *1 condizioney

logical *1 si

logical *1 cluster

parameter ( ipoints = 100 000 )

parameter (maxelements =20 000) ! dimensione vettori vari

dimension kubo ( 0:mmm,0:mmm )

dimension cluster ( 0:mmm,0:mmm )

dimension icx ( 1:maxelements )

dimension icy ( 1:maxelements )

dimension ielements ( 1:maxelements )

dimension inewx ( 1:maxelements )

dimension inewy ( 1:maxelements )

dimension inewxsur ( 1:maxelements )

dimension inewysur ( 1:maxelements )

dimension isx ( 1:ipoints )

dimension isy ( 1:ipoints )

ifine = 0

do jx = 0 , mmm

do jy = 0 , mmm

if ( kubo (jx ,jy) ) then

ifine = ifine + 1

isx (ifine) = jx

isy (ifine) = jy

153

5.3. PERCOLAZIONE CAPITOLO 5. SISTEMI MONTE CARLO

Figura 5.5: Cluster di connessione

154

CAPITOLO 5. SISTEMI MONTE CARLO 5.3. PERCOLAZIONE

end if

end do

end do

if ( ifine.gt.ipoints) then

stop ’aumentare ipoints !’

end if

icluster = 0

1015 continue

icluster = icluster +1

do jj=1,ifine

ix = isx (jj)

iy = isy (jj)

if (kubo(ix,iy)) then

icx (1) = ix

icy (1) = iy

kubo(ix,iy ) =.false.

inewx (1 ) = ix

inewy (1 ) = iy

go to 102

end if

end do

icluster = icluster -1

go to 1016

102 continue

mmm1 = mmm + 1

number = 1

new = 1

103 continue

iadd = 0

cccccccccccccccc caricamento primi vettori di nuovi punti ccccccccccccc

do j=1,new

do jcasi = 1, 4

if ( jcasi.eq. 1 ) then

isux = 0

isuy = 1

elseif ( jcasi.eq. 2 ) then

isux = +1

isuy = 0

elseif ( jcasi.eq. 3 ) then

155

5.3. PERCOLAZIONE CAPITOLO 5. SISTEMI MONTE CARLO

isux = 0

isuy = -1

elseif ( jcasi.eq. 4 ) then

isux = -1

isuy = 0

end if

newx = inewx ( j ) + isux

newy = inewy ( j ) + isuy

if (( newx.eq.-1).or.(newx.eq.mmm1)) then

go to 1050

end if

if (( newy.eq.-1).or.(newy.eq.mmm1)) then

go to 1050

end if

if (( isux.eq.0).and.(isuy.eq.0 )) then

go to 1050

end if

c x boundaries problem

if (kubo(newx,newy)) then

kubo(newx,newy)=.false.

iadd = iadd +1

inewxsur ( iadd ) = newx

inewysur ( iadd ) = newy

end if

1050 continue

end do

end do

c determinzaione 4 vicini ccccccccccccccc

if ( iadd.eq.0 ) then

ielements ( icluster) = number

cc inseriamo il loop sul calcolo della percolazione

cc scanning su y

condizionea = .false.

condizioneb = .false.

do j= 1 , number

if (icy(j).eq. 0 ) then

condizionea = . true .

end if

if (icy(j).eq. mmm ) then

156

CAPITOLO 5. SISTEMI MONTE CARLO 5.3. PERCOLAZIONE

condizioneb = . true .

end if

end do

condizioney=.false.

if ((condizionea).and.(condizioneb)) then

condizioney =.true.

end if

cc scanning su x

condizionea = .false.

condizioneb = .false.

do j= 1 , number

if (icx(j).eq. 0 ) then

condizionea = . true .

end if

if (icx(j).eq. mmm ) then

condizioneb = . true .

end if

end do

condizionex =.false.

if ((condizionea).and.(condizioneb)) then

condizionex =.true.

if (condizionex.and.condizioney) then

si =.true.

do jjj = 1 , number

iixx = icx ( jjj )

jjyy = icy ( jjj )

cluster ( iixx , jjyy ) =.true.

end do

end if

end if

go to 1015

end if

do j = 1,iadd

number = number + 1

icx (number ) = inewxsur (j)

icy (number ) = inewysur (j)

end do

do j = 1,iadd

inewx ( j ) = inewxsur (j)

157

5.4. CATENE DI VORONOI CAPITOLO 5. SISTEMI MONTE CARLO

inewy ( j ) = inewysur (j)

end do

new = iadd

go to 103

ccccccccc end great loop cccccccccccccccccccccccccccccccccccccc

1016 continue

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

type *,’clusters veri=’,icluster

return

end

5.4 Catene di Voronoi

Questi diagrammi prendono il nome da un matematico russo che nel1908 fece un lavoro ( ovviamente analitico ) sui poligoni irregolari cheriempiono il piano. Il problema e poi stato ripreso da svariati auto-ri ognuno con in mente una applicazione di diverso tipo: strutturadei cristalli, struttura della cellula, raffreddamento e conseguente fram-mentazione delle rocce, struttura dell’ universo,etc. Il problema in duedimensioni consiste nel trovare una volta che e dato un insieme di punti,detti nuclei, una partizione dello spazio tale che ogni nucleo abbia unsuo dominio e le frontiere, ovverosia le zone equidistanti fra un nucleoe l’ altro, siano marcate.

Per riassumere chiamiamo diagrammi di Voronoi il set dei punti taliche le distanze fra i due nuclei piu vicini siano uguali. Elenchiamo unaserie possibile di disposizioni di nuclei nel piano classificabile in basealla distribuzione di probalita di partenza oppure in base all’ algoritmogenerante:

1. generazione random su x ed y,

2. generazione gaussiana polare su r e θ,

3. numeri pseudo-random o di Sobol,

4. autovalori di matrice complessa,

5. punti frattali.

158

CAPITOLO 5. SISTEMI MONTE CARLO 5.4. CATENE DI VORONOI

Figura 5.6: Diagramma di Voronoi con nuclei random su x e y

159

5.4. CATENE DI VORONOI CAPITOLO 5. SISTEMI MONTE CARLO

Figura 5.7: Diagramma di Voronoi con nuclei gaussiani

160

CAPITOLO 5. SISTEMI MONTE CARLO 5.4. CATENE DI VORONOI

Figura 5.8: La testa della margherita

161

5.4. CATENE DI VORONOI CAPITOLO 5. SISTEMI MONTE CARLO

Si tratta poi di trovare tramite uno scanning degli N ×N punti dellattice quelli equidistanti dai due nuclei piu vicini per quanto riguarda ilati e dai tre nuclei piu vicini per quanto riguarda i vertici. Riportiamonelle due figure 5.6 e 5.7 un tipico esempio di tessellazione quando inuclei sono distribuiti in maniera random su x y ed in maniera gaussia-na . Le applicazioni pratiche sono svariate e come esempio riportiamola struttura logaritmica della testa dei fiori tipo margherita o girasole.Per generare i nuclei usiamo la formula di Vogel:

Φ = n× 137.5 r = c√n , (5.3)

dove

1. n e il numero di punti in senso crescente dal centro. Questo el’ inverso dell’ eta dell’ organismo vegetale in questione.

2. Φ e l’ angolo di referenza in gradi del punto n in coordinate polari.

3. r e la distanza fra il centro ed il punto n.

Il risultante diagramma di Voronoi viene riportato nella figura 5.8Per maggiori dettagli sui diagrammi di Voronoi consultare il libro [OBS 92] .L’ implementazione pratica dello scanning avviene attraverso la SUB-ROUTINE MATRICE.

subroutine

, matrice (vec,ipunti,side,matrix,nodes,mll,muu,mmm,traspo)

implicit real * 8 (a-h,o-z)

logical * 1 matrix,yes,nodes,lato1,lato2

parameter (nnnn=1000)

dimension matrix (0:mmm,0:mmm )

dimension nodes (0:mmm,0:mmm )

dimension traspo (1:20 )

dimension vec (1:ipunti,1:2)

dimension ads0 (1:nnnn)

dimension ads (1:nnnn)

dimension awork (1:nnnn)

dimension index (1:nnnn)

dimension indexwork (1:nnnn)

dimension nodesx (1:nnnn)

dimension nodesy (1:nnnn)

162

CAPITOLO 5. SISTEMI MONTE CARLO 5.4. CATENE DI VORONOI

delta = side/mmm

delta0 = side/mll

rad2 = dsqrt (2.0 D+00)

correz = 1.

adelta = delta * correz

adelta0 = delta0 * correz

deltamezzi = delta /2.0

delta0mezzi= delta0/2.0

do jm=0,mmm

do jn=0,mmm

matrix (jm,jn)=.false.

nodes (jm,jn)=.false.

end do

end do

numero = 0

do jy0 =0 , mll -1

do jx0 =0 , mll -1

xx = jx0 * delta0 + delta0mezzi

yy = jy0 * delta0 + delta0mezzi

do jj =1,ipunti

ax = xx - vec(jj,1)

ay = yy - vec(jj,2)

ads0 (jj) = dsqrt ( ax*ax + ay*ay )

end do

ifail =0

call m01caf ( ads0,1,ipunti,’A’,ifail) ! riordino del vettore in

! senso ascendente

adiff0 = dabs (ads0(2)- ads0 (1))

if (adiff0.le.( adelta0* rad2)) then

do jy =jy0*muu , jy0*muu + muu-1

do jx =jx0*muu , jx0*muu + muu-1

xx = jx * delta + deltamezzi

yy = jy * delta + deltamezzi

do jj =1,ipunti

ax =xx - vec(jj,1)

ay =yy - vec(jj,2)

ads (jj) = dsqrt ( ax*ax + ay*ay )

end do

ifail =0

163

5.4. CATENE DI VORONOI CAPITOLO 5. SISTEMI MONTE CARLO

call m01ajf (ads,awork,index,indexwork,ipunti,ipunti,ifail)

adiff = dabs (ads(2)- ads (1))

cccccccccccc beginning check on all the square vertex ccccccccccccccccc

if (adiff.le.(adelta*rad2)) then

x1 = vec(index(1),1)

y1 = vec(index(1),2)

x2 = vec(index(2),1)

y2 = vec(index(2),2)

xc = xx

yc = yy

yes=.false.

call retta90 (x1,y1,x2,y2,xc,yc,deltamezzi,yes)

if (yes) then

numero = numero +1

matrix (jx,jy)=.true.

cccccccccc beginning nodes cccccccccccccccccccccccccccccccccccccccc

lato1 =.false.

lato2 =.false.

x1 = vec(index(1),1)

y1 = vec(index(1),2)

x2 = vec(index(3),1)

y2 = vec(index(3),2)

xc = xx

yc = yy

yes=.false.

call retta90 (x1,y1,x2,y2,xc,yc,deltamezzi,yes)

if (yes) then

lato1=.true.

end if

x1 = vec(index(2),1)

y1 = vec(index(2),2)

x2 = vec(index(3),1)

y2 = vec(index(3),2)

xc = xx

yc = yy

yes=.false.

call retta90 (x1,y1,x2,y2,xc,yc,deltamezzi,yes)

if (yes) then

lato2=.true.

164

CAPITOLO 5. SISTEMI MONTE CARLO 5.4. CATENE DI VORONOI

end if

if ((lato1).and.(lato2)) then

nodes (jx,jy)=.true.

endif

ccccccccccccccccccccccccc end nodes ccccccccccccccccccccccccccc

end if

end if

end do

end do

end if

end do

end do

traspo (1) = adelta

traspo (2) = numero

return

end

Questa subroutine fa uso della SUBROUTINE RETTA90.

subroutine retta90 (x1,y1,x2,y2,xc,yc,deltamezzi,yes)

implicit real *8 (a-h,o-z)

logical * 1 yes,uno,due

yes=.false.

uno=.false.

due=.false.

deltay = y2 -y1

deltax = x2 -x1

am = deltay/deltax

100 continue

c determinazione primo coefficiente angolare

amprimo = -1.0/am

c determinazione secondo coefficiente angolare a 90 gradi

ymezzi = (y1+y2)/2.0

xmezzi = (x1+x2)/2.0

165

5.4. CATENE DI VORONOI CAPITOLO 5. SISTEMI MONTE CARLO

x = xc + deltamezzi

y = ymezzi + amprimo * ( x - xmezzi)

c determinazione prima retta

ylower = yc - deltamezzi

yupper = yc + deltamezzi

if ((y.gt.ylower).and.(y.le.yupper)) then

uno=.true.

end if

x = xc - deltamezzi

y = ymezzi + amprimo * ( x - xmezzi)

c determinazione prima retta

if ((y.gt.ylower).and.(y.le.yupper)) then

uno=.true.

end if

y = yc + deltamezzi

x = (y-ymezzi)/amprimo + xmezzi

xlower = xc - deltamezzi

xupper = xc + deltamezzi

if ((x.gt.xlower).and.(x.le.xupper)) then

due=.true.

end if

y = yc - deltamezzi

x = (y-ymezzi)/amprimo + xmezzi

if ((x.gt.xlower).and.(x.le.xupper)) then

due=.true.

end if

if ((uno).or.(due)) then

yes = .true.

end if

return

end

166

Appendice A

Funzioni intrinseche

Riportiamo il nome generico delle funzioni intrinseche attualmente inuso sui compilatori di FORTRAN 77 insieme ad una breve spiegazionedel loro significato.

ABS (numero) Una funzione che ritorna il valore assoluto dell’ ar-gomento. Il valore assoluto di un numero complesso, (X,Y), e il suovalore reale (X**2 + Y**2)**(1/2).

ACOS (numero)Una funzione che ritorna l’ arcocoseno dell’ argo-mento in radianti. Il valore assoluto dell’ argomento deve essere minoreo uguale a 1.

ACOSD ( numero reale)Una funzione che ritorna l’ arcocosenodell’ argomento in gradi. Il valore dell’ argomento deve essere compresofra 0 (escluso) e 1 (incluso).

AIMAG (numero complesso)Una funzione che ritorna la parteimmaginaria di un numero complesso. L’ argomento deve essere di tipoCOMPLEX*8. Il risultato e di tipo REAL*4.

AINT ( numero reale)Una funzione che ritorna l’ intero piu gran-de il cui valore assoluto non supera il valore assoluto dell’ argomento eha lo stesso segno dell’ argomento.

AMAX0 (numero, numero,...)Una funzione che ritorna il mag-giore dei valori specificato nella lista degli argomenti.

AMIN0 (numero, numero,...) Una funzione che ritorna il mi-nore dei valori specificati nella lista degli argomenti.

ANINT ( numero reale) Una funzione che ritorna il valore dell’ in-tero piu vicino al valore dell’ argomento.

167

APPENDICE A. FUNZIONI INTRINSECHE

ASIN ( numero reale) Una funzione che ritorna l’ arcoseno dell’ ar-gomento in radianti. Il valore assoluto dell’ argomento deve essereminore o uguale a 1.

ASIND ( numero reale) Una funzione che ritorna l’ arcosenodell’ argomento in gradi. Il valore dell’ argomento deve essere compresofra 0 (escluso) e 1 (incluso).

ATAN ( numero reale) Una funzione che ritorna l’ arcotangentedell’ argomento in radianti. Il valore assoluto dell’ argomento deveessere maggiore di 0.

ATAND ( numero reale) Una funzione che ritorna l’ arcotan-gente dell’ argomento in gradi. Il valore assoluto dell’ argomento deveessere maggiore di 0.

ATAN2 ( numero reale, numero reale) Una funzione che ri-torna l’ arcotangente del quoziente dei due argomenti in radianti. Seambedue gli argomenti sono zero, il risultato e indefinito. Se il primoargomento e positivo, il risultato e positivo. Se il primo argomento e ne-gativo, il risultato e negativo. Se il primo argomento e zero, il risultatoe zero. Se il secondo argomento e zero, il valore assoluto del’ risultatoe pi/2. Il range del risultato e compreso tra -pi e pi.

ATAN2D ( numero reale, numero reale) Una funzione cheritorna l’ arcotangente del quoziente dei due argomenti in gradi. Seambedue gli argomenti sono zero, il risultato e indefinito. Se il primoargomento e positivo, il risultato e positivo. Se il primo argomentoe negativo, il risultato e negativo. Se il primo argomento e zero, ilrisultato e zero. Se il secondo argomento e zero, il valore assolutodel’ risultato e 90. Il range del risultato e compreso tra -180 e 180gradi.

BTEST (intero, posizione) Una funzione che ritorna il valorelogico.true. se il bit con l’ intero specificato dalla posizione e settato a1. I Il bit di ordine basso e posizionato a 0.

CHAR (intero) Una funzione che ritorna il carattere associatocon il valore ASCII specificato dall’ argomento. L’ argomento deve es-sere BYTE, LOGICAL*1, LOGICAL*2, LOGICAL*4, INTEGER*2,o INTEGER*4 come tipo di dato. Il risultato deve essere di tipoCHARACTER. Il valore di input deve essere nel range da 0 a 255.

CMPLX (numero [,numero]) Una funzione che converte uno odue argomenti nel tipo COMPLEX*8. Se viene specificato un argo-mento, questo e convertito nella parte reale del valore del complesso e

168

APPENDICE A. FUNZIONI INTRINSECHE

la parte immaginaria diventa zero. Se sono specificati due argomenti, ilprimo e convertito nella parte reale del valore del complesso e il secondonella parte immaginaria del valore del complesso. I due argomenti chesono specificati devono essere dello stesso tipo.

CONJG (numero complesso) Una funzione che ritorna il com-plesso coniugato dell’ argomento. Se l’ argomento e (X,Y), il suocomplesso coniugato e (X,-Y).

COS (numero) Una funzione che ritorna il coseno dell’ argomento.L’ argomento deve essere in radianti.

COSD (numero) Una funzione che ritorna il coseno dell’ argo-mento. L’ argomento deve essere in gradi.

COSH ( numero reale) Una funzione che ritorna il coseno iper-bolico dell’ argomento.

DATE (data) Una subroutine che mette la data corrente nel suoargomento. L’ argomento a 9-byte puo essere una stringa carattere,una variable, vettore, o un elemento di vettore di ogni tipo. La data eun valore ASCII formattato come segue: gg-mmm-aa.

DBLE (numero) Una funzione che converte l’ argomento in unvalore REAL*8.

DCMPLX (numero [,numero]) Una funzione che converte unoo due argomenti nel tipo COMPLEX*16. Se viene specificato un argo-mento, questo e convertito nella parte reale del valore del complesso ela parte immaginaria diventa zero. Se sono specificati due argomenti,il primo e convertito nella parte reale del valore del complesso e il se-condo argomento nella parte immaginaria del valore del complesso. Idue argomenti che sono specificati devono essere dello stesso tipo.

DFLOAT (integer) Una funzione che converte l’ argomento in unvalore REAL*8.

DIM (numero, numero) Una funzione che ritorna il valore delprimo argomento meno il minimo (MIN) dei due argomenti.

DIMAG (immaginario) Una funzione che ritorna la parte im-maginaria di un numero complesso. L’ argomento deve essere di tipoCOMPLEX*16. Il risultato e del tipo REAL* 8.

DPROD (numeroreal*4, numeroreal*4) Una funzione che ri-torna il prodotto di due valori REAL*4 come un valore REAL*8.

DREAL (numero complesso) Una funzione che ritorna la par-te reale di un numero complesso. L’ argomento deve essere di tipoCOMPLEX*16. Il risultato e di tipo REAL*8.

169

APPENDICE A. FUNZIONI INTRINSECHE

EXP (esponente) Una funzione che ritorna e**X, dove X e ilvalore dell’ argomento.

FLOAT (integer) Una funzione che converte l’ argomento in unvalore REAL*4.

IABS (numero) Una funzione che ritorna il valore assoluto di unargomento intero. Il valore assoluto di un numero complesso, (X,Y), eil valore reale (X**2 + Y**2)**(1/2).

IAND (intero, intero) Una funzione che produce il LOGICALdegli argomenti su una base del tipo bit per bit.

IBCLR (intero, posizione) Una funzione che ritorna il valore delprimo argomento con il bit specificato settato a 0. Il bit di ordine bassoe in posizione 0.

IBITS (intero, posizione di partenza, lunghezza ) Una fun-zione che ritorna il valore del bit specificato dalla posizione di partenzae dalla posizione finale. Il bit di ordine piu basso e in posizione 0.

IBSET (intero, posizione) Una funzione che ritorna il valore delprimo argomento con il bit specificato settato a 0. Il bit di ordine piubasso e in posizione 0.

ICHAR (carattere) Una funzione che ritorna il valore ASCIIdell’ argomento. L’ argomento deve essere una espressione characterdi lunghezza uno. Il risultato e un tipo di dato INTEGER*4.

IDATE ( mese, giorno, anno ) Una subroutine che ritorna itre valori rappresentanti la data corrente. Gli argomenti devono esseredefiniti come interi o vettore di interi. Il mese e rappresentato con ilnumero del mese (1 -12). Il giorno e rappresentato come il giorno delmese. L’ anno e rappresentato dalle ultime due cifre dell’ anno.

IDIM (numero, numero) Una funzione che ritorna il valore delprimo argomento meno il minimo (MIN) dei due argomenti.

IDINT (numero) Una funzione che ritorna l’ intero piu grande ilcui valore assoluto non eccede il valore assoluto dell’ argomento e ha lostesso segno dell’ argomento. L’ argomento deve essere di tipo real * 8.

IDNINT ( numero reale) Una funzione che ritorna il valoredell’ intero piu vicino al valore dell’ argomento. L’ argomento deveessere di tipo real * 8.

IEOR (intero, intero) Una funzione che riporta un OR esclusivosu argomenti di tipo bit per bit.

IFIX (real*4) Una funzione che converte un numero reale in unintero

170

APPENDICE A. FUNZIONI INTRINSECHE

INDEX (stringa, sub-stringa ) Una funzione che cerca una strin-ga per la sua prima comparsa di una sottostringa e ritorna la posizionedi partenza della sottostringa come un valore INTEGER*4.

INT (numero) Una funzione che ritorna il piu grande intero ilcui valore assoluto non eccede il valore assoluto dell’ argomento e halo stesso segno dell’ argomento. La funzione ritorna un valore INTE-GER*4 se il commando / I4 e in funzione; altrimenti ritorna un valoreINTEGER*2.

IOR (intero, intero) Una funzione che produce un OR logicodell’ argomento sulla base di ragionamenti bit per bit.

IQINT (numero) Una funzione che ritorna l’ intero piu grande ilcui valore assoluto non eccede il valore assoluto dell’ argomento ed halo stesso segno dell’ argomento. L’ argomento deve essere di tipo real* 16.

ISIGN (valore, segno )Una funzione che assegna il segno delsecondo argomento al valore assoluto del primo.

LEN (carattere)Una funzione che ritorna il numero di caratteridell’ argomento. L’ argomento deve essere un’ espressione CHARAC-TER. Il risultato e un valore INTEGER*4.

LGE (carattere, carattere )Una funzione che ritorna il valo-re.true. se la prima stringa carattere e piu lunga o uguale alla secondastringa carattere. Gli argomenti devono essere del tipo CHARACTERmentre il risultato e un valore LOGICAL*4.

LGT (carattere, carattere )Una funzione che ritorna il valo-re.true. se la prima stringa carattere e maggiore della seconda stringacarattere. Gli argomenti devono essere del tipo CHARACTER mentreil risultato e un valore LOGICAL*4.

LLE (carattere, carattere )Una funzione che ritorna il valo-re.true. se la prima stringa carattere e minore o uguale alla secondastringa carattere. Gli argomenti devono essere del tipo CHARACTERmentre il risultato e un valore LOGICAL*4.

LLT (carattere, carattere )Una funzione che ritorna il valo-re.true. se la prima stringa carattere e minore alla seconda stringacarattere. Gli argomenti devono essere del tipo CHARACTER mentreil risultato e un valore LOGICAL*4.

LOG (numero) Una funzione che ritorna il logaritmo naturale(base e) dell’ argomento. L’ argomento deve essere maggiore di zero.

171

APPENDICE A. FUNZIONI INTRINSECHE

LOG10 (numero) Una funzione che ritorna il logaritmo comune(base 10 ) dell’ argomento. L’ argomento deve essere maggiore di zero.

MAX (numero, numero,...) Una funzione che ritorna il maggioredei valori specificati nella lista degli argomenti.

MIN (numero, numero,...) Una funzione che ritorna il minoredei valori specificati nella lista degli argomenti.

MOD (dividendo, divisore) Una funzione che divide il primoargomento per il secondo e ritorna il resto.

NINT ( numero reale) Una funzione che ritorna il valore dell’ in-tero piu vicino al valore dell’ argomento.

QEXT (numero)Una funzione che converte l’ argomento in unvalore REAL*16.

QFLOAT (intero) Una funzione che converte un valore intero inun REAL*16.

RAN (seme ) Una funzione che ritorna un valore un numeroREAL*4 fra 0.0 (incluso) e 1.0 (escluso). L’ argomento deve essereuna variabile o elemento di vettore INTEGER*4. Per ottenere i risul-tati migliori inizializzate con un valore grande e dispari per la primavolta. Per generare sets differenti di valori random inizializzate il semeogni volta con un valore diverso.

REAL (numero) Una funzione che converte l’ argomento in unvalore real.

SIZEOF ( argomento ) La funzione ritorna il numero di bytes dimemoria usata dall’ argomento.

SECNDS ( numero reale) Una subroutine che ritorna il numerodi secondi da mezzanotte meno il valore dell’ argomento. L’ argomentodeve essere di tipo REAL*4. Il valore di ritorno e del tipo REAL*4. Iltempo di ritorno e accurato a 1/100 di secondo.

SIGN (valore, segno ) Una funzione che assegna il segno delsecondo argomento al valore assoluto del primo.

SIN (numero) Una funzione che ritorna il seno dell’ argomento.L’ argomento deve essere in radianti.

SIND (numero) Una funzione che ritorna il seno dell’ argomento.L’ argomento deve essere in gradi.

SINH (numero) Una funzione che ritorna il seno iperbolico dell’ ar-gomento che deve essere espresso in radianti.

SINHD (numero) Una funzione che ritorna il seno iperbolicodell’ argomento che deve essere espresso in gradi.

172

APPENDICE A. FUNZIONI INTRINSECHE

SQRT (numero) Una funzione che ritorna la radice quadratadell’ argomento. Il risultato e il valore principale, con la parte realemaggiore o uguale a zero. Se la parte reale e zero, il risultato e il valoreprincipale, con l’ immaginario maggiore o uguale a zero.

TAN ( numero reale) Una funzione che ritorna la tangente dell’ ar-gomento. L’ argomento deve essere in radianti.

TAND ( numero reale) Una funzione che ritorna la tangentedell’ argomento. L’ argomento deve essere in gradi.

TANH ( numero reale) Una funzione che ritorna la tangente iper-bolica dell’ argomento. L’ argomento deve essere espresso in radianti.

TANHD ( numero reale) Una funzione che ritorna la tangenteiperbolica dell’ argomento. L’ argomento deve essere espresso in gradi.

TIME (tempo) Una subroutine che ritorna il tempo attuale informato ASCII 24-h. L’ argomento deve essere una variabile, vettori,elemento di vettore o sottostringa di 8 bytes. Il tempo e formattatocome segue: hh:mm:ss.

ZEXT (intero) Una funzione che ritorna il valore dell’ argomento,zero esteso.

173

APPENDICE A. FUNZIONI INTRINSECHE

174

Appendice B

Subroutine di interessecomune

Riportiamo in questa appendice i listati delle subroutines adoperate piufrequentemente nei programmi precedenti.

B.1 CLEANTOP

Questa subroutine pulisce lo schermo e porta il cursore a sinistra in altosui terminali.

subroutine cleantop

implicit real*4 (a-h,o-z)

c

c this program clears the screen , then changes the screen from

c dark to light

print *,char (27),’[2J’

c clear the screen

print *,char (27),’[0;0H’

c bring the cursor at the top left

c write (6,*)char(27),’[?5h’

c change the screen background to light

return

end

175

B.2. MAXMIN APPENDICE B. SUBROUTINE DI INTERESSE COMUNE

B.2 MAXMIN

Questa subroutine ritorna il massimo ( XXMAX) e il minimo ( XXMIN) di un vettore XX di dimensioni IVEC.

subroutine maxmin (xx,ivec,xxmax,xxmin)

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c c

c this subroutine is producing a the max and the minimum of the c

c vector xx

c ivec :== vector dimension c

c c

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

implicit real * 4 (a-h,o-z)

dimension xx(1:ivec)

xxmax = xx (1)

xxmin = xx (1)

do j = 1, ivec

if (xx(j).gt.xxmax) then

xxmax = xx(j)

end if

if (xx(j).lt.xxmin) then

xxmin = xx(j)

end if

end do

return

end

176

Appendice C

Sistema operativo UNIX

Si entra su una macchina UNIX tramite la comunicazione dilogin ePassword. Compiuta questa prima operazione comparira un promptche dipende dal sistema operativo usato. Se nel file .login inserite ilcomando set prompt = ” % ” comparira il segno di percentuale adaccompagnarvi nelle vostre avventure. Riportiamo adesso i comandiprincipali che interessano a noi.

C.1 Alcuni comandi UNIX

% man nomecomando (C.1)

E molto importante per chi non ha tempo o voglia di leggersi i manualiin quanto permette di vedere cosa corrisponde alla voce nomecomando.

% man − k stringa (C.2)

Serve a cercare una particolare stringa inserita in un comando. Laschermata successiva visualizza il comando coinvolto e la posizionecontestuale della stringa.

% ls (C.3)

Questo comando riporta su schermo tutti i files.

% ls − l (C.4)

Per vedere la loro grandezza ( numero di bytes ) e la data in cui furonocreati.

177

C.2. EDITOR APPENDICE C. SISTEMA OPERATIVO UNIX

% ls − l − a (C.5)

Serve a vedere i file nascosti (quelli che incominciano con punto tipo.netscape , .csrhc, .login etc....).

% cat miofile (C.6)

Permette un evince del file indicato sul terminale adoperato.

% rm miofile (C.7)

Cancella il file indicato; ricordiamo che in questo sistema operativoesiste solo l’ ultima versione del file creato.

% mv vecchiofile nuovofile (C.8)

Cambia il nome del file indicato in quello nuovo.

% cp miofile altrofile (C.9)

Copia il file indicato in uno nuovo.

% exit (C.10)

Serve a chiudere la sessione interattiva.

< Ctrl/c > (C.11)

Serve ad interrompere il processo interattivo in corso e ritornare alsistema operativo.

C.2 Editor

E il mezzo fondamentale per avere accesso ai files nuovi oppure a quelliesistenti e ne abbiamo svariati. Ci limitiamo a quelli piu semplici chesono stati anche testati

% joe miofile (C.12)

Si puo accedere all’ help on line dell’ editor con il comando ( in editing) < Ctrl +K +H >

% pico miofile (C.13)

Questo e l’ editor che adotta la finestra del noto sistema di postaelettronica che risponde al nome di PINE.

178

APPENDICE C. SISTEMA OPERATIVO UNIXC.3. COMPILAZIONE+LINK/ESECUZIONE

C.3 Compilazione+link/esecuzione

Per compilare ed eseguire un programma, ad esempio miofile.f, si dannoi seguenti due comandi:

% f77 − o miofile.exe miofile.f (C.14)

% ./miofile.exe (C.15)

Che sono in ordine comando per la compilazione + link ed esecuzione.

Se invece vogliamo solamente controllare che il programma noncontenga eventuali errori di compilazione, diamo il comando

% f77 − g − c miofile.f . (C.16)

C.4 Shell Script

Spesso puo essere utile abbreviare una intera serie di comandi trami-te un solo file dettoscript. Come esempio riportiamo uno script percompilare ed eseguire un programma che chiameremocalcola.

calcola#!/bin/shclearecho ”Inserire nome del file”read nomeecho ”ho inserito”echo $nomerm $nome.exeecho ”la data e l’ ora sono:‘date‘ ”f77 -o $nome.exe $nome.for./$nome.exe

Volendo adesso eseguire questo file di comandi dovremo prima cam-biare le protezioni con il comandochmod u+x calcola e poi il coman-do ./calcola. Cercate tramite il manuale di capire le varie operazionicompiute da questa sequenza.

179

C.5. LE LIBRERIE APPENDICE C. SISTEMA OPERATIVO UNIX

C.5 Le librerie

Se volete crearvi una vostra libreria personale dovete prima trasformarei file sorgente in file oggetto e poi collezionarli in un unica libreria.Riportiamo uno script che genera una libreria chiamata mia lib.f77 -c *.far rcv libmia lib.a *.oranlib libmia lib.arm *.o

Controllate sulle pagine del manuale il significato delle varie opera-zioni. Se avete il file della libreria potete linkare ad essa il programmamain con il seguente comando:f77 main.f -L. -lmia lib

180

Appendice D

I caratteri ASCII

I sistemi collegati con i computers memorizzano i caratteri sotto seriedi bits, usualmente lunghi 7 bits oppure 8. Riportiamo di seguito lasequenza dei 128 caratteri ASCII ( American Standard Code for In-formation Interchange ) in notazione decimale. Ricordiamo che esistepure la rappresentazione ottale e quella esadecimale.

181

APPENDICE D. I CARATTERI ASCII

Tabella D.1: Tabella caratteri ASCII

NUL 0 DLE 16 SP 32 0 48 @ 64 P 80 ‘ 96 p 112SOH 1 DC1 17 ! 33 1 49 A 65 Q 81 a 97 q 113STX 2 DC2 18 ” 34 2 50 B 66 R 82 b 98 r 114ETX 3 DC3 19 # 35 3 51 C 67 S 83 c 99 s 115EOT 4 DC4 20 $ 36 4 52 D 68 T 84 d 100 t 116ENQ 5 NAK 21 % 37 5 53 E 69 U 85 e 101 u 117ACK 6 SY N 22 & 38 6 54 F 70 V 86 f 102 v 118BEL 7 ETB 23 ′ 39 7 55 G 71 W 87 g 103 w 119BS 8 CAN 24 ( 40 8 56 H 72 X 88 h 104 x 120HT 9 EM 25 ) 41 9 57 I 73 Y 89 i 105 y 121LF 10 SUB 26 ∗ 42 : 58 J 74 Z 90 j 106 z 122V T 11 ESC 27 + 43 ; 59 K 75 [ 91 k 107 { 123FF 12 FS 28 , 44 ¡ 60 L 76 \ 92 l 108 | 124CR 13 GS 29 − 45 = 61 M 77 ] 93 m 109 } 125SO 14 RS 30 . 46 > 62 N 78 ˆ 94 n 110 ˜ 126SI 15 US 31 47 ? 63 O 79 95 o 111 DEL 127

182

Bibliografia

[Ba 88] M.Barnsley,“Fractals everywhere ”,prima edizione, Academic Press ,Inc. , Boston(1988).

[Do 83] R. Doretti,R. Farabone,“Dal FORTRAN IV al FORTRAN 77”,prima edizione, Gruppo Editoriale Jackson , Milano (1983).

[El 85] T.M.R. Ellis,“Programmazione strutturata in FORTRAN 77 ”,prima edizione, Zanichelli , Bologna (1985).

[Ol 92] D. Oliver,“Fractal Visions ”,prima edizione, SAMS PUBLISHING , Carmel (1992).

[OBS 92] Okabe,A. ,Boots,B.,Sugihara,K. ,“Spatial Tessellations Concepts and Applications of Voronoi Dia-grams”,prima edizione, John Wiley & Sons , Chichester (1992).

[Or 86] A. Orsi Palamara ,“Programmare in FORTRAN 77 ”,prima edizione, Levrotto e Bella , Torino (1986).

[PS88] H.O. Peitgen,S. Saupe,“The Science of Fractal Image ”,prima edizione, Springer Verlag , Londra (1988).

[Pa 83] C. Page,“FORTRAN 77 ”,prima edizione, Pitman Publishing LTD , Londra (1983).

[REC 92] Press W.H.,Teukolsky S.A., Vetterling W.T.,Flannery B.P.,“Numerical Recipes in Fortran”,seconda edizione, Cambridge University Press, Cambridge (1992).

183

BIBLIOGRAFIA BIBLIOGRAFIA

[PL 90] P. Prusinkiewicz, A. Lindenmayer, “The Algorithmic Beautyof Plants”,prima edizione, Springer Verlag, New York (1990).

[TF 92] Stauffer D., Aharony A., “Introduction to Percolation Theo-ry”,seconda edizione, Taylor & Francis , London (1992).

[Vi 92] T.vicsek, “Fractal Growth Phenomena”,seconda edizione, Word scientific, Singapore (1992).

184

Indice analitico

AACCEPT istruzione, 40aggregazione-ballistica, 142aggregazione-DLA, 133algoritmi-Montecarlo, 133argomenti-subroutine, 20ASCII-caratteri, 179assegnazioni, 27ASSIGN istruzione, 28

BBACKSPACE istruzione, 41

CCALL istruzione, 19caratteri, 2catene-Voronoi, 156CHARACTER variabile, 7, 8, 11CHARACTER istruzione, 5CLOSE istruzione, 38colonne , 3commento, 4COMMON istruzione, 24COMPLEX variabile, 7COMPLEX istruzione, 5concatenazione, 11continuazione, 4CONTINUE istruzione, 33

DDATA istruzione, 27

DIMENSION istruzione, 24dimensione-vettori, 20distribuzione

esponenziale, 98gaussiana, 99Student, 95

DO istruzione, 31DO-WHILE istruzione, 33DOUBLE PRECISION variabile,

7DOUBLE PRECISION istruzione,

5

EELSE istruzione, 30ELSEIF istruzione, 30END DO istruzione, 32ENDFILE istruzione, 41END istruzione, 5ENTRY istruzione, 21EQUIVALENCE istruzione, 25espressioni, 9Etichette, 5EXTERNAL istruzione, 25

Ffile, 36fit

retta, 105format

/, 44

185

INDICE ANALITICO INDICE ANALITICO

Aw, 44BN, 44BZ, 44Dw.d , 43Ew.d, 43Fw.d, 42Gw.d , 43IW, 42KP, 44Lw, 43nHc1..cn , 43SP, 44SS, 44stringa, 44Tln, 44Tn, 44Trn, 44

FORMAT istruzione, 42FORTRAN, 1FUNCTION istruzione, 16funzione

beta incompleta, 93coefficiente-binomiale, 87errori, 91fattoriale, 86gamma, 88gamma-incompleta, 88generatori-interi, 97generatori-sistema, 95numeri-random, 95p-chiquadro, 92

funzione intrinseca, 165ABS, 165ACOS, 165ACOSD, 165AIMAG, 165AINT, 165AMAX0, 165AMIN0, 165

ANINT, 165ASIN, 166ASIND, 166ATAN, 166ATAN2, 166ATAN2D, 166ATAND, 166BTEST, 166CHAR, 166CMPLX, 167CONJG, 167COS, 167COSD, 167COSH, 167DATE, 167DBLE, 167DCMPLX, 167DFLOAT, 167DIM, 167DIMAG, 167DPROD, 167DREAL, 167EXP, 168FLOAT, 168IABS, 168IAND, 168IBCLR, 168IBITS, 168IBSET, 168ICHAR, 168IDATE, 168IDIM, 168IDINT, 168IDNINT, 168IEOR, 168IFIX, 168INDEX, 169INT, 169IOR, 169

186

INDICE ANALITICO INDICE ANALITICO

IQINT, 169ISIGN, 169LEN, 169LGE, 169LGT, 169LLE, 169LLT, 169LOG, 169LOG10, 170MAX, 170MIN, 170MOD, 170NINT, 170QEXT, 170QFLOAT, 170RAN, 170REAL, 170SECNDS, 170SIGN, 170SIN, 170SIND, 170SINH, 170SINHD, 170SIZEOF, 170SQRT, 171TAN, 171TAND, 171TANH, 171TANHD, 171TIME, 171ZEXT, 171

Gglobali, 4GO TO istruzione, 29

IIF istruzione, 30IMPLICIT istruzione, 22

IMPLICIT istruzione, 4INQUIRE istruzione, 39INTEGER istruzione, 5INTEGER variabile, 6INTRINSIC istruzione, 26

Llocali, 4LOGICAL variabile, 7LOGICAL istruzione, 5

Mmain, 14mapping

Feigenbaum, 112Henon, 109Julia, 119regolabile, 120

memoria, 6momenti-campione, 101

asimmetria, 102curtosi, 102deviazione-assoluta, 101deviazione-standard-empirica, 101valor-medio, 101varianza-empirica, 101

OOPEN istruzione, 37

PPARAMETER istruzione, 22parole-chiave, 5PAUSE istruzione, 34Percolazione, 148PRINT istruzione, 40priorita-istruzioni, 15PROGRAM, 16

187

INDICE ANALITICO INDICE ANALITICO

PROGRAMMACAVEN, 60GAUSS, 74PARABO, 50PENDO, 52RETTA, 47testc2, 80testt, 79

RREAD istruzione, 39REAL istruzione, 5REAL variabile, 6record, 35relazioni, 12relazioni-logiche, 12RETURN istruzione, 21REWIND istruzione, 40

SSAVE istruzione, 26scheda, 3sottostringhe, 8STOP istruzione, 34SUBROUTINE

AVEVAR, 105BA, 142BETA, 88BETACF, 94BETAI, 93BICO, 87CLEANTOP, 173DUETUTTO, 119ERF, 92ERFC, 92EXPDEV, 99FACTLN, 87FACTRL, 86FAI-CLUSTER, 156

FAI-PUNTI, 151FEIGEN, 115FELCI, 126FIT, 106G05DYFZAN, 98GAMMA, 85GAMMLN, 85GAMMP, 89GAMMQ, 89GASDEV, 100GCF, 90GSER, 90GX, 132HENON10, 112MAKEPIANTA3, 130MATRICE, 160MAXMIN, 174MOMENT, 102RAN0, 97RETTA90, 163SCRIVIDATI, 41TTTEST, 104UNO, 118WO, 134

SUBROUTINE istruzione, 18

Ttest

medie, 103tipi, 5TYPE istruzione, 40

UUNIX

cp, 176exit, 176joe, 176ls, 175man, 175

188

INDICE ANALITICO INDICE ANALITICO

mv, 176rm, 176

Vvettori, 8

WWRITE istruzione, 39

189