mockus: 'dati dell'idrogramma unitario di Mockus
DATA 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4,4.1,4.2,4.3,4.4,4.5,4.6,4.7,4.8,4.9,5:
't/ta DATA 0.030,0.100,0.190,0.310,0.470,0.660,0.820,0.930,0.990,1.000,0.990,0.930,0.860,0.780,0.680,0.560,0.460,0.390,0.330,0.280,0.244,0.207,0.177,0.147,0.127,0.107,0.092,0.077,0.066,0.055,0.048,0.040,0.035,0.029,0.025,0.021,0.018,0.015,0.013,0.011,0.009,0.008,0.007,0.006,0.006,0.005,0.004,0.003,0.001,0.000:
' t/tp
coefficienti:
DATA 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,20,22,24,26,28,30,32:
'N volte tc!. DATA 0.1,0.1,0.1,0.2,0.3,0.2,0.3,0.3,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,1,1,1,1,1,1,1,1:
'discretizzazione dei coeff. di Mockus.
m
AS INTEGER '[numero naturale] conta il numero di interazioni fino al raggiungimento di "d". mdt
AS SINGLE '[ore] conta la progessione delle ore fino al raggiungimento di "d". tminuti
AS SINGLE '[minuti] trasforma mdt in minuti. h1
AS SINGLE '[mm] valore h della curva di possibilità climatica. i1
AS SINGLE '[mm/ora] intensità di pioggia. i2
AS SINGLE '[mm/ora] intensità di pioggia secondo lo ietogramma Chicago o costante. DH
AS SINGLE '[mm] pioggia in ogni istante temporale dt. H2
AS SINGLE '[mm] pioggia cumulata ad ogni istante temporale. he
AS SINGLE '[mm] piggia efficace cumulata ad ogni istante temporale. Dhe
AS SINGLE '[mm] pioggia efficace ad ogni istante temporale. qmSUqp
AS SINGLE '[-] il valore di qm/qp corrispondente a tm/ta e' = al valore di q/qp che corrisponde a t/ta=tm/ta.
DIM SHARED S2!
'[mm] contenuto idrico massimo del terreno. DIM SHARED qp!
'[mc/s] portata al colmo dell'idrogramma unitario. DIM SHARED k!
'[-] coeff. della curva di possibilità climatica tado dal sito della Regione. DIM SHARED a2!
'[mm/d^n] coeff. della curva di possibilità climatica. DIM SHARED n1!
'[-] coeff. della curva di possibilità climatica. DIM A1!
'[mq] superficie bacino idrografico. DIM L!
'[m] lunghezza asta principale. DIM HmaxL!
'[m] punto più alto dell'astra principale. DIM HminL!
'[m] punto più basso dell'asta principale. DIM s1!
'[-] pendenza media del bacino non espressa in percentuale ma come 0,x. DIM CNII!
'[-] tabellato. DIM CNIII
'[-] discende da CNII. DIM tl!
'[ore] tempo di ritardo.
A1! = 2641902
k! = 2.34
a2! = 14.6653
n1! = 0.45596
L! = 3300
HmaxL! = 1981
HminL! = 1062
s1! = 0.415
CNII! = 70
'INPUT A1!
'INPUT k!
'INPUT a2!
'INPUT n1!
'INPUT L!
'INPUT HmaxL!
'INPUT HminL!
'INPUT s1!
'INPUT CNII!
'INPUT d1!
CNIII! = (23 * CNII!) / (10 + 0.13 * CNII!)
tl! = 0.342 * ((L! / 1000) ^ 0.8 / (100 * s1!) ^ 0.5) * (1000 / CNIII! - 9) ^ 0.7 'formula di Mockus.
S2! = 25.4 * (1000 / CNIII! - 10)
Ia! = 0.1 * S2! 'coeff.=0.03-0.2.
tc! = tl! / 0.6
ta! = tl! / 0.9
qp! = 0.208 * ((A1! / 1000000) / ta!)
IF i%
<= 24 THEN READ coefficienti
(i%
).d: coefficienti
(i%
).d
= coefficienti
(i%
).d
* tc!
'PRINT "Ietogramma Chicago/costante: 1/2> "
'DO
' _LIMIT 60 'limit
'keypress$ = INKEY$
'LOOP UNTIL keypress$ = "1" OR keypress$ = "2"
'visualizzaieto% = VAL(keypress$)
visualizzaieto% = 1
CALL grafici
(massimi1
(visualizzaieto%
, 24, 1).ore
, massimi2
(visualizzaieto%
).ore
, massimi2
(visualizzaieto%
).portata
) 'valori passati per disegnare il grafico. 'CALL risultati
'-----------------------------------------matrice1---------------------------------------------------------------------------------------------------------------------------------
z% = 0
z% = z% + 1
dt!(ieto%, z%) = coefficienti(z%).tSUta * ta!
m%(ieto%, z%) = coefficienti(z%).d / (coefficienti(z%).tSUta * ta!) 'determina il numero di passi temporali necessari ad arrivare allla durata d considerata con il paso temporale dt!.
matrice1(ieto%, z%, i%, 1).MOCKUStSUta = mockus(i%).MOCKUStSUta
matrice1(ieto%, z%, i%, 1).MOCKUSqSUqp = mockus(i%).MOCKUSqSUqp
matrice1(ieto%, z%, i%, 1).m = i%
matrice1(ieto%, z%, i%, 1).mdt = matrice1(ieto%, z%, i%, 1).m * dt!(ieto%, z%)
matrice1(ieto%, z%, i%, 1).tminuti = matrice1(ieto%, z%, i%, 1).mdt * 60
SELECT CASE (_ROUND(10 * (i%
* coefficienti
(z%
).tSUta
))) / 10 'crea il vettore dell'idrogramma unitario di Mockus con i multipli del coeff. z% considerato (0.1,0.2,0.3,0.5,1). CASE IS <= 5 'valore massimo di t/ta di Mockus matrice1
(ieto%
, z%
, i%
, 1).tmSUta
= (_ROUND(10 * (i%
* coefficienti
(z%
).tSUta
))) / 10 '=matrice1(i%).m * matrice1(i%).mdt / ta! matrice1(ieto%, z%, i%, 1).tmSUta = 0
n% = 1
DO 'crea il vettore dell'idrogramma unitario Mockus prendendo i valori corrispondenti di ai coefficienti. IF matrice1
(ieto%
, z%
, i%
, 1).tmSUta
= matrice1
(ieto%
, z%
, n%
, 1).MOCKUStSUta
THEN matrice1(ieto%, z%, i%, 1).qmSUqp = matrice1(ieto%, z%, n%, 1).MOCKUSqSUqp
n% = n% + 1
matrice1(ieto%, z%, i%, 1).qm = matrice1(ieto%, z%, i%, 1).qmSUqp * qp!
FOR i%
= 1 TO m%
(ieto%
, z%
) matrice1(ieto%, z%, i%, 1).h1 = k! * a2! * matrice1(ieto%, z%, i%, 1).mdt ^ n1!
matrice1(ieto%, z%, i%, 1).i1 = matrice1(ieto%, z%, i%, 1).h1 / dt!(ieto%, z%)
matrice1(ieto%, z%, i%, 1).i1 = (matrice1(ieto%, z%, i%, 1).h1 - matrice1(ieto%, z%, i% - 1, 1).h1) / dt!(ieto%, z%)
FOR i%
= 1 TO m%
(ieto%
, z%
) matrice1(ieto%, z%, i%, 1).i2 = matrice1(ieto%, z%, m%(ieto%, z%), 1).h1 / matrice1(ieto%, z%, m%(ieto%, z%), 1).mdt
FOR i%
= 1 TO m%
(ieto%
, z%
) matrice1(ieto%, z%, i%, 1).DH = matrice1(ieto%, z%, i%, 1).i2 * dt!(ieto%, z%)
matrice1(ieto%, z%, i%, 1).H2 = matrice1(ieto%, z%, i%, 1).DH
matrice1(ieto%, z%, i%, 1).H2 = matrice1(ieto%, z%, i%, 1).DH + matrice1(ieto%, z%, i% - 1, 1).H2
matrice1(ieto%, z%, i%, 1).he = 0
matrice1(ieto%, z%, i%, 1).he = (matrice1(ieto%, z%, i%, 1).H2 - Ia!) ^ 2 / (matrice1(ieto%, z%, i%, 1).H2 - Ia! + S2!)
matrice1(ieto%, z%, i%, 1).Dhe = matrice1(ieto%, z%, i%, 1).he
matrice1(ieto%, z%, i%, 1).Dhe = matrice1(ieto%, z%, i%, 1).he - matrice1(ieto%, z%, i% - 1, 1).he
'-------------------------------matrice2-------------------------------------------------------------------------------------------------------------------------
FOR n%
= 1 TO 50 'colonna matrice2!(ieto%, z%, i%, n%) = matrice1(ieto%, z%, i% - n% + 1, 1).Dhe * matrice1(ieto%, z%, n%, 1).qm
'-------------------------------idrogramma------------------------------------------------------------------------------------------------------------------
FOR n%
= 1 TO 50 'colonna idrogramma(ieto%, z%, i%, 1).portata = idrogramma(ieto%, z%, i%, 1).portata + matrice2!(ieto%, z%, i%, n%)
idrogramma(ieto%, z%, i%, 1).ore = matrice1(ieto%, z%, i%, 1).mdt
idrogramma(ieto%, z%, i%, 1).minuti = matrice1(ieto%, z%, i%, 1).tminuti
'-------------------------------Massimi-----------------------------------------------------
i% = 1
'massimi1(ieto%, z%, 1).portata = idrogramma(ieto%, z%, 1, 1).portata
IF massimi1
(ieto%
, z%
, 1).portata
< idrogramma
(ieto%
, z%
, i%
, 1).portata
THEN massimi1(ieto%, z%, 1).portata = idrogramma(ieto%, z%, i%, 1).portata
massimi1(ieto%, z%, 1).ore = idrogramma(ieto%, z%, i%, 1).ore
massimi1(ieto%, z%, 1).minuti = idrogramma(ieto%, z%, i%, 1).minuti
i% = i% + 1
' IF massimi2(ieto%).ore < massimi1(visualizzaieto%, z%, 1).ore THEN massimi2(ieto%).ore = massimi1(visualizzaieto%, z%, 1).ore
IF massimi2
(ieto%
).portata
< massimi1
(visualizzaieto%
, z%
, 1).portata
THEN massimi2(ieto%).portata = massimi1(visualizzaieto%, z%, 1).portata
massimi2(ieto%).ore = massimi1(visualizzaieto%, z%, 1).ore
'---------------------------------------------------------------------------------------------------------------------------------------------------
'Questa subroutine crea gli ietogrammi chicago per ogni iterazione z%.
'Ad ogni nuovo dt! pone l'ultimo valore di pioggia come ultimo, il penultimo come primo, il terzultimo come penultimo, il quartultimo come secondo e così via.
'Quando vi è un valore di z% che ha un dt! non nuovo, cerca la precedente iterazione z% con lo stesso dt! e copia i suoi valori dello ietogramma chicago in quello dell'iterazione z% in corso. I valori rimanenti sono riempiti
'procedendo come prima.
'Quindi per esempio l'iterazione z%=1 corrisponde a d=1tc con passo temporale dt!. L'iterazione z%=2 corrisponde a d=2tc con il medesimo passo temporale. Per cui lo ietogramma chicago per d=2tc e' dato dallo stesso ietogramma
'chicago che si ha per d=1tc, con l'aggiunta, ai suoi estremi, dei nuovi valori. Il meccanismo si ripete per z%=3, che corrisponde a d=3tc, con il medesimo dt!. Per z%=4, corrispondente a d=4tc, invece, il passo temporale dt!
'cambia.
'Quanto detto e' necessario per evitare piccole discrepanze che potrebbero comportare il fatto che alcuni massimi risultino localmente poco piu' bassi del precedente, anche quando il trend generale e' chiaramente di aumento dei
'massimi all'aumentare della durata della pioggia. In questo modo si fa l'ipotesi che all'aumentare della durata della pioggia, discretizzata con uno stesso passo temporale dt!, ogni ietogramma contenga a sua volta il precedente.
'In tal modo, se si condiera per esempio d=2tc, si assume che l'altezza di pioggia corrispondente a 1tc, segua lo ietogramma relativo a d=1tc, cosi' che l'altezza di pioggia rimamente sia eclusivamente dovuta ai nuovi valori.
n% = 1
IF z%
= 1 THEN test%
= 0:
EXIT DO 'permette il calcolo dello ietogramma chicago della prima iterazione. test%=0 IF dt!
(ieto%
, z%
) = dt!
(ieto%
, z%
- n%
) THEN ' se trova un'iterazione z% precedente con lo stesso dt! di quella in corso, esegue i comandi successivi test% = 1
k% = m%(ieto%, z%)
i% = 0
DO 'inserisce i nuovi valori dello ietogramma chicago dell'iterazione z% in corso. matrice1(ieto%, z%, m%(ieto%, z%) - i%, 1).i2 = matrice1(ieto%, z%, k%, 1).i1
k% = k% - 1
IF k%
= m%
(ieto%
, z%
- n%
) THEN EXIT DO 'se ha completato l'inserimento dei nuovi valori, passa oltre. matrice1(ieto%, z%, i% + 1, 1).i2 = matrice1(ieto%, z%, k%, 1).i1
k% = k% - 1
i% = i% + 1
IF matrice1
(ieto%
, z%
, i%
, 1).i2
= 0 THEN i%
= i%
- 1 k% = 1
DO 'inserisce lo ietogramma chicago dell'iterazione z% precedente con lo stesso dt!, nello ietogramma dell'iterazione z% in corso. i% = i% + 1
matrice1(ieto%, z%, i%, 1).i2 = matrice1(ieto%, z% - n%, k%, 1).i2
k% = k% + 1
EXIT DO 'siccome ha completato lo ietogramma dell'iterazione z% in corso, non serve incrementare n%, quindi passa oltre. test%=1. ELSE ' se non trova l'iterazione z% precedente con lo stesso dt! di quella in corso, incrementa n% e ripete la ricerca. test% = 0
n% = n% + 1
LOOP UNTIL n%
= z%
- 1 'se non ha trovato iterazioni z% precedenti con lo stesso dt!, allora passa oltre e test%=0. CASE IS = 1 'esce dalla subroutine, in quanto lo ietogramma e' definito. CASE IS = 0 'calcola ex-novo lo ietogramma relativo all'iterazione z% in corso in quanto non vi e' un'iterazioni z% precedente con lo stesso dt! dell'iterazione z% in corso. k% = m%(ieto%, z%)
i% = 0
matrice1(ieto%, z%, m%(ieto%, z%) - i%, 1).i2 = matrice1(ieto%, z%, k%, 1).i1
k% = k% - 1
matrice1(ieto%, z%, i% + 1, 1).i2 = matrice1(ieto%, z%, k%, 1).i1
k% = k% - 1
i% = i% + 1
'-----------------------------------------------------------------------------------------------------------------------------------------------------------------
SUB grafici
(ore!
, oreM!
, QM!
) 'ore!: ore corrispondenti al massimo maggiore.
'oreM!: ore corrispondenti alla portata massima.
'QM!: portata massima.
DIM L%
, H%
'dimensioni del grafico.
DIM dx%
, dy%
'origine degli assi dentro SCREEN. DIM dx!
, dy!
'posizione dell'origine degli assi dentro WINDOW. DIM x%
, y%
'conversione nelle coordinate di SCREEN delle coordinate dentro WINDOW corrispondenti a determinati valori.
H%
= _DESKTOPHEIGHT * 2 / 3 ' imposta l'altezza dello SCREEN pari ai 2/3 dello schermo visualizzabile, cosi' che su ogni PC sia visualizzabile. L% = H% * 1.62 'imposta la larghezza dello SCREEN in modo che L%xH% sia un rettangolo aureo.
dx% = 39
dy% = H% - 1 - dx%
taccaX% = ore! \ 5 ': PRINT taccaX%: SLEEP
taccaY% = QM! \ 5
' SCREEN grafico&
LINE (0, 0)-(L%
- 1, H%
- 1), R
, B
WINDOW (0, 0)-(ore!
* 1.15, QM!
* 1.15) 'scala ascisse e ordinate in base alle ore e alla portata, estendendole del 15%. LINE (dx!
, QM!
* 1.1)-(dx!
, dy!
):
LINE -(ore!
* 1.1, dy!
), bianco&
'disegna gli assi in modo che si estendano del 10% oltre valori massimi, ma senza arrivare alla massima estensione delle coordinate impostate da WINDOW: cosi' si ha che gli assi non toccano la cornice rossa e i valori massimi che non arrivano al limite degli assi. PSET (dx!
, QM!
* 1.1!
), bianco&:
DRAW "F20":
PSET (dx!
, QM!
* 1.1), bianco&:
DRAW "G20" PSET (ore!
* 1.1, dy!
), bianco&:
DRAW "G20":
PSET (ore!
* 1.1, dy!
), bianco&:
DRAW "H20" colore& = R&
colore& = G&
colore& = B&
CIRCLE (dx!
+ massimi1
(visualizzaieto%
, z%
, 1).ore
, dy!
+ massimi1
(visualizzaieto%
, z%
, 1).portata
), PMAP(3, 2), R&
PAINT (dx!
+ massimi1
(visualizzaieto%
, z%
, 1).ore
+ PMAP(0.5, 2), dy!
+ massimi1
(visualizzaieto%
, z%
, 1).portata
+ PMAP(0.5, 2)), R&
' _LIMIT limit
LINE (dx!
, dy!
)-(dx!
+ idrogramma
(visualizzaieto%
, z%
, i%
, 1).ore
, dy!
+ idrogramma
(visualizzaieto%
, z%
, i%
, 1).portata
), colore&
LINE (dx!
+ idrogramma
(visualizzaieto%
, z%
, i%
- 1, 1).ore
, dy!
+ idrogramma
(visualizzaieto%
, z%
, i%
- 1, 1).portata
)-(dx!
+ idrogramma
(visualizzaieto%
, z%
, i%
, 1).ore
, dy!
+ idrogramma
(visualizzaieto%
, z%
, i%
, 1).portata
), colore&
' _PUTIMAGE ((_DESKTOPWIDTH - L%) \ 2, (_DESKTOPHEIGHT - H%) \ 2), grafico&, schermo&
LINE (dx!
, dy!
)-(dx!
+ massimi1
(visualizzaieto%
, z%
, 1).ore
, dy!
+ massimi1
(visualizzaieto%
, z%
, 1).portata
), giallo&
LINE (dx!
+ massimi1
(visualizzaieto%
, z%
- 1, 1).ore
, dy!
+ massimi1
(visualizzaieto%
, z%
- 1, 1).portata
)-(dx!
+ massimi1
(visualizzaieto%
, z%
, 1).ore
, dy!
+ massimi1
(visualizzaieto%
, z%
, 1).portata
), giallo&
LINE (dx!
, dy!
+ QM!
)-(dx!
+ oreM!
, dy!
+ QM!
), bianco&
, , 61440 LINE -(dx!
+ oreM!
, dy!
), bianco&
, , 61440 x%
= PMAP(dx!
+ oreM!
, 0) i% = 1
WHILE i%
* taccaX%
<= ore!
LINE (dx!
+ i%
* taccaX%
, PMAP(dy%
+ 5, 3))-(dx!
+ i%
* taccaX%
, PMAP(dy%
- 5, 3)) x%
= PMAP(dx!
+ i%
* taccaX%
, 0) i% = i% + 1
i% = 1
WHILE i%
* taccaY%
<= QM!
LINE (PMAP(dx%
- 5, 2), dy!
+ i%
* taccaY%
)-(PMAP(dx%
+ 5, 2), dy!
+ i%
* taccaY%
) y%
= PMAP(dy!
+ i%
* taccaY%
, 1) i% = i% + 1
' DO
' CLS 2
' PRINT "Digitare le ore e i minuti (<= "; _TRIM$(STR$(INT(ore!))); " ore e "; _TRIM$(STR$(INT((ore! - INT(ore!)) * 60))); " minuti) di cui si vuole conoscere la portata di picco."
'INPUT "Ore: ", ore2!
'INPUT "Minuti: ", minuti%
'ore2! = ore2! + minuti% / 60
'LOOP WHILE ore2! > ore!
ore2! = 3
x%
= PMAP(dx!
+ ore2!
, 0) SCREEN grafico&
'------------------------------------------------------------------------------------>????????????????????????????????????? y% = y% - 1
' _LIMIT limit
LINE (dx!
+ ore2!
, dy!
)-(dx!
+ ore2!
, PMAP(y%
, 3)), bianco&
, , 65280 ' _FREEIMAGE grafico&
'------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
PRINT "dt="; dt!
(visualizzaieto%
, z%
);
";"; dt!
(visualizzaieto%
, z%
) * 60 PRINT "visualizzaieto%="; visualizzaieto%
PRINT "m%="; m%
(visualizzaieto%
, z%
) PRINT "N="; coefficienti
(z%
).d
/ tc!
PRINT "n="; coefficienti
(z%
).tSUta
' _LIMIT limit
PRINT _TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).m
)), _TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).mdt
)), _TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).tminuti
)), _TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).h1
)),_
_TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).i1
)),_TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).i2
)),_TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).DH
)), _TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).H2
)),_
_TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).he
)),_TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).Dhe
)),_TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).MOCKUStSUta
)),_
_TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).MOCKUSqSUqp
)),_TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).tmSUta
)), _TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).qmSUqp
)),_TRIM$(STR$(matrice1
(visualizzaieto%
,z%
,i%
, 1).qm
)) ' _LIMIT limit
PRINT matrice2!
(visualizzaieto%
, z%
, i%
, n%
) ' _LIMIT limit
'PRINT _TRIM$(STR$(idrogramma(visualizzaieto%, n%, i%, 1).ore))
'PRINT _TRIM$(STR$(idrogramma(visualizzaieto%, n%, i%, 1).minuti))
LOCATE i%
, (n%
- 16) * 15 - 14 'PRINT _TRIM$(STR$(idrogramma(visualizzaieto%, n%, i%, 1).ore))
'PRINT _TRIM$(STR$(idrogramma(visualizzaieto%, n%, i%, 1).minuti))
' _LIMIT limit
PRINT massimi1
(visualizzaieto%
, z%
, 1).ore
, massimi1
(visualizzaieto%
, z%
, 1).minuti
, massimi1
(visualizzaieto%
, z%
, 1).portata
PRINT massimi2
(visualizzaieto%
).ore; massimi2
(visualizzaieto%
).portata