-
Notifications
You must be signed in to change notification settings - Fork 44
/
Copy pathdlev.gs
134 lines (109 loc) · 2.76 KB
/
dlev.gs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
*
* Help is in the end of this script
*
function dlev( args )
_version='0.01r2'
if( args = '' )
help()
return
endif
***** Default value *****
var = ''
dvdlev = ''
loglev = 0
***** Arguement *****
i = 1
while(1)
arg = subwrd( args, i )
i = i + 1;
if( arg = '' ); break; endif
while(1)
*** option
if( arg = '-loglev' ); loglev=1; break; endif
*** int, max, min
if( var = '' ); var=arg; break; endif
if( var != '' & dvdlev = '' ); dvdlev=arg; break; endif
say 'syntax error : 'arg
return
endwhile
endwhile
'q dims'
line = sublin( result, 4 )
flag = subwrd( line, 3 )
if( flag != 'varying' ); return; endif
zmin = subwrd( line, 11 )
zmax = subwrd( line, 13 )
zminpp = zmin + 1
zmaxmm = zmax - 1
* Half level
'set z 'zmin' 'zmaxmm
if( loglev = 0 )
'dvp = ( 'var'(z+1) - 'var' ) / ( lev(z+1) - lev )'
'levp = 0.5 * ( lev + lev(z+1) )'
else
'dvp = ( 'var'(z+1) - 'var' ) / ( log(lev(z+1)) - log(lev) )'
'levp = 0.5 * ( log(lev) + log(lev(z+1)) )'
endif
* 'dv = 0 * 'var
* z = zmin
* while( z <= zmaxmm )
* 'set z 'z+1
* 'varpp = 'var
* 'set z 'zmin' 'zmaxmm
** 'dv = dv + maskout(varpp,z)'
* z = z + 1
* endwhile
*exit
'set z 'zminpp' 'zmax
if( loglev = 0 )
'dvm = ( 'var' - 'var'(z-1) ) / ( lev - lev(z-1) )'
'levm = 0.5 * ( lev(z-1) + lev )'
else
'dvm = ( 'var' - 'var'(z-1) ) / ( log(lev) - log(lev(z-1)) )'
'levm = 0.5 * ( log(lev(z-1)) + log(lev) )'
endif
* Full level
'set z 'zminpp' 'zmaxmm
if( loglev = 0 )
'w = ( lev - levm ) / ( levp - levm )'
'out = w * dvp + (1-w) * dvm'
else
'w = ( log(lev) - levm ) / ( levp - levm )'
'out = w * dvp + (1-w) * dvm'
'out = 'out' / lev'
endif
* Output
'set z 'zmin' 'zmax
if( dvdlev = '' )
'd out'
else
dvdlev' = out'
endif
return
*
* help
*
function help()
say ' Name:'
say ' dlev '_version' - differentiate variable with respect to lev'
say ' '
say ' Usage:'
say ' dlev var [var_out] [-loglev]'
say ''
say ' var : variable to be differentiated'
say ' var_out : result variable'
say ' if var_out is not specified,'
say ' result is displayed on the screen instead'
say ' -loglev : use log(lev) instead of lev'
say ''
say ' Note:'
say ' [arg-name] : specify if needed'
say ' (arg1 | arg2) : arg1 or arg2 must be specified'
say ' Without -loglev, var_out = d(var)/d(lev)'
say ' With -loglev, var_out = d(var)/dlog(lev) / lev'
say ' Half-level value is calculated for accuracy (2nd-order)'
say ''
say ' Copyright (C) 2009 Chihiro Kodama'
say ' Distributed under GNU GPL (http://www.gnu.org/licenses/gpl.html)'
say ''
return