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 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
|
subroutine isoval(a,lda,m,n,c,path,npath,maxp,ierr,iw,job)
c! but
c Etant donnee une matrice A telle que A(l,k)=f(X(l),Y(k)) ou
c f est une fonction de R X R dans R, ce sous programme recherche
c les lignes d'isovaleurs (relatives a la valeur c) de la
c tabulation reguliere ( X(l+1)-X(l)=DX,Y(k+1)-Y(k)=DY )
c de f donnee par A
c! liste d'appel
c subroutine isoval(a,lda,m,n,c,path,npath,maxp,ierr,iw,job)
c
c double precision a(lda,n),c,path(2,maxp)
c integer iw(m*n-n)
c integer lda,m,n,maxp,ierr,job
c
c a : tableau contenant la tabulation de f
c lda : nombre de lignes du tableau a
c m : nombre de lignes effectif de a
c n : nombre de colonnes de a
c c : valeur pour la quelle on cherche les isovaleurs
c path : contient en sortie la description des isovaleurs:
c path=[path ,...., path ]
c 1 npath
c ou pathi a la structure suivante :
c [np x1,...xnp]
c [0 y1,...ynp]
c si :
c np est la longueur de l'isovaleur
c xj,yj les coordonnees interpolees des points de l'isovaleur
c npath : le nombre de courbes disjointes d'isovaleurs
c maxp : la dimension maximale admise pour le tableau path
c ierr : indicateur d'erreur
c 0 : ok
c 1 : nombre de points decrivant les isovaleur > maxp
c iw : tableau de travail
c job : flag permettant d'indiquer au programme si la fonction f
c est definie sur l'ensemble des points de la matrice A
c job=0 : f definie partout
c job=1 : f n'est pas definie au points A(i,j) si
c iw(i+(j-1)m)=0
c! origine
c programme par S Steer a partir de la macro scilab de Carey Bunks
c 1990
c
c Corrige par C Bunks pour isovaleurs qui sont exactement egales
c a des entrees de la matrice 6 mars 1991.
c!
c Copyright INRIA
double precision a(*),c,path(2,*)
integer iw(*)
integer lda,m,n,maxp,ierr,job
c
logical pend
integer north,south,east,west,dir
data north/0/,south/1/,east/2/,west/3/
c
ierr=0
c
kv=1
kh=kv+m*n
kw=kh+m*(n-1)
c
c perturb values which are exactly equal to the level value
do 5 ip=1,n*m
if(a(ip).eq.c) a(ip)=a(ip)+1d-14
5 continue
c
c make horizontal and vertical edge matrices for level value
if(job.eq.0) then
ih=kh-1
ia=-lda
do 11 k=1,n-1
ia=ia+lda
do 10 l=1,m
ih=ih+1
iw(ih)=1
if((a(ia+lda+l)-c)*(a(ia+l)-c).lt.0.0d0) iw(ih)=-1
10 continue
11 continue
c
iv=kv-1
ia=-lda
do 13 k=1,n
ia=ia+lda
do 12 l=1,m-1
iv=iv+1
iw(iv)=1
if((a(ia+l+1)-c)*(a(ia+l)-c).lt.0.0d0) iw(iv)=-1
12 continue
13 continue
else
c
kj=kv
c
ih=kh-1
ia=-lda
ij=kj-1
do 21 k=1,n-1
ia=ia+lda
do 20 l=1,m
ij=ij+1
ih=ih+1
iw(ih)=1
if(iw(ij)*iw(ij+m).eq.0) goto 20
if((a(ia+lda+l)-c)*(a(ia+l)-c).lt.0.0d0) iw(ih)=-1
20 continue
21 continue
c
iv=kv-1
ia=-lda
ij=kj-1
do 23 k=1,n
ia=ia+lda
do 22 l=1,m-1
iv=iv+1
ij=ij+1
iw(iv)=1
if(iw(ij)*iw(ij+1).eq.0) goto 22
if((a(ia+l+1)-c)*(a(ia+l)-c).lt.0.0d0) iw(iv)=-1
22 continue
ij=ij+1
23 continue
endif
c
npath=0
kpath0=1
kpath=0
c
c search pathes (starting with boundaries)
c
c horizontal boundaries
c northern border
ih=kh-m
ia=1-lda
do 31 ic=1,n-1
ih=ih+m
ia=ia+lda
if(iw(ih).lt.0) then
kpath=1
path(2,kpath0+1)=1.0d0
path(1,kpath0+1)=ic+(c-a(ia))/(a(ia+lda)-a(ia))
i=1
j=ic
dir=north
pend=.false.
30 call isova0(a,lda,m,n,path(1,kpath0+1),kpath,i,j,dir,pend,
$ iw(kh),iw(kv),c)
if(kpath0+kpath.ge.maxp) goto 999
if(.not.pend) goto 30
if(kpath.gt.1) then
path(1,kpath0)=kpath
path(2,kpath0)=0
kpath0=kpath0+1+kpath
npath=npath+1
endif
kpath=0
endif
31 continue
c
c southern border
c
ih=kh+(m-1)+(n-1)*m
ia=1+(m-1)+(n-1)*lda
do 41 ic=n-1,1,-1
ih=ih-m
ia=ia-lda
if(iw(ih).lt.0) then
kpath=1
path(2,kpath0+1)=m
path(1,kpath0+1)=ic+(c-a(ia))/(a(ia+lda)-a(ia))
i=m-1
j=ic
dir=south
pend=.false.
40 call isova0(a,lda,m,n,path(1,kpath0+1),kpath,i,j,dir,pend,
$ iw(kh),iw(kv),c)
if(kpath0+kpath.ge.maxp) goto 999
if(.not.pend) goto 40
if(kpath.gt.1) then
path(1,kpath0)=kpath
path(2,kpath0)=0
kpath0=kpath0+kpath+1
npath=npath+1
endif
kpath=0
endif
41 continue
c
c vertical boundaries
c
c eastern border
iv=kv-1+(n-1)*(m-1)
ia=(n-1)*lda
do 51 ir=1,m-1
iv=iv+1
ia=ia+1
if(iw(iv).lt.0) then
kpath=1
path(2,kpath0+1)=ir+(c-a(ia))/(a(ia+1)-a(ia))
path(1,kpath0+1)=n
i=ir
j=n-1
dir=east
pend=.false.
50 call isova0(a,lda,m,n,path(1,kpath0+1),kpath,i,j,dir,pend,
$ iw(kh),iw(kv),c)
if(kpath0+kpath.ge.maxp) goto 999
if(.not.pend) goto 50
if(kpath.gt.1) then
path(1,kpath0)=kpath
path(2,kpath0)=0
kpath0=kpath0+1+kpath
npath=npath+1
endif
kpath=0
endif
51 continue
c
c western border
c
iv=kv+m-1
ia=m
do 61 ir=m-1,1,-1
iv=iv-1
ia=ia-1
if(iw(iv).lt.0) then
kpath=1
path(2,kpath0+1)=ir+(c-a(ia))/(a(ia+1)-a(ia))
path(1,kpath0+1)=1.0d0
i=ir
j=1
dir=west
pend=.false.
60 call isova0(a,lda,m,n,path(1,kpath0+1),kpath,i,j,dir,pend,
$ iw(kh),iw(kv),c)
if(kpath0+kpath.ge.maxp) goto 999
if(.not.pend) goto 60
if(kpath.gt.1) then
path(1,kpath0)=kpath
path(2,kpath0)=0
kpath0=kpath0+kpath+1
npath=npath+1
endif
kpath=0
endif
61 continue
c
c all the rest
c
ih=kh-1
ia1=1-lda
do 72 ic=1,n-1
ia1=ia1+lda
ia=ia1
ih=ih+1
do 71 ir=2,m-1
ih=ih+1
ia=ia+1
if(iw(ih).lt.0) then
kpath=1
path(2,kpath0+1)=ir
path(1,kpath0+1)=ic+(c-a(ia))/(a(ia+lda)-a(ia))
i=ir
j=ic
dir=north
pend=.false.
70 call isova0(a,lda,m,n,path(1,kpath0+1),kpath,i,j,dir,
$ pend,iw(kh),iw(kv),c)
if(kpath0+kpath.ge.maxp) goto 999
if(.not.pend) goto 70
if(kpath.gt.1) then
path(1,kpath0)=kpath
path(2,kpath0)=0
kpath0=kpath0+kpath+1
npath=npath+1
endif
kpath=0
endif
71 continue
ih=ih+1
72 continue
return
999 ierr=1
end
|