1 : C INTERP SOURCE GOUNAND 12/09/11 21:15:03 7496
2 : SUBROUTINE INTERP
3 : C
4 : C=======================================================================
5 : C
6 : C Opérateur IPOL
7 : C
8 : C----------------------------
9 : C Phrases d'appel (GIBIANE) :
10 : C----------------------------
11 : C
12 : C OBJET2 = 'IPOL' OBJET1 | LISTREE1 LISTREE2 | ;
13 : C | EVOL1 |
14 : C
15 : C OBJET3 = 'IPOL' TAB1 FLOT1 ;
16 : C
17 : C
18 : C OBJET2 = 'IPOL' OBJET1 NUA1 | ('GAUSS') | ;
19 : C | 'RATIO' |
20 : C
21 : C------------------------
22 : C Opérandes et résultat :
23 : C------------------------
24 : C
25 : C Syntaxe 1
26 : C
27 : C LISTREE1 : liste des abscisses T de la fonction
28 : C
29 : C LISTREE2 : liste des ordonnées F(T) de la fonction
30 : C
31 : C EVOL1 : Evolution représentant la fonction
32 : C
33 : C OBJET1 : objet contenant les valeurs à interpoler sur F
34 : C (type FLOTTANT, CHPOINT, LISTREEL, MCHAML)
35 : C
36 : C OBJET2 : valeurs de la fonction interpolée aux points donnés
37 : C (même type que OBJET1)
38 : C
39 : C Syntaxe 2
40 : C
41 : C TAB1 : table de sous-type 'RESULTAT'
42 : C
43 : C FLOT1 : valeur à interpoler sur F
44 : C
45 : C OBJET3 : valeurs de la fonction interpolée au point donné
46 : C (type CHPOINT ou MCHAML)
47 : C
48 : C Syntaxe 3
49 : C
50 : C NUA1 : objet nuage représentant la fonction à interpoler
51 : C
52 : C=======================================================================
53 : C
54 : C Remarques
55 : C
56 : C Les listes LISTREE1 et LISTREE2 doivent se correspondre
57 : C
58 : C L'évolution EVOL1 doit être élémentaire
59 : C
60 : C ATTENTION : la liste des abscisses donnéee est supposée monotone
61 : C
62 : C=======================================================================
63 : C
64 : C Auteur :
65 : C Création :
66 : C Modifications :
67 : C PM le 03/10/2006 : lecture d'une évolution en entrée
68 : C 19/03/2007 : vérification que les listes ont une longueur > 1
69 : C SG 22/10/2009 : option SPLINE (biblio : numerical recipes)
70 : C=======================================================================
71 :
72 : IMPLICIT INTEGER(I-N)
73 : IMPLICIT REAL*8(A-H,O-Z)
74 : CHARACTER*8 TYPE1,TYPE2,TYPEI,TYPO1,TYPO2,TYPOBJ
75 : CHARACTER*8 CHARIN,CHARRE
76 : LOGICAL LOGIN,LOGRE
77 : CHARACTER*4 CLE(1)
78 : DATA CLE/'SPLI'/
79 : CHARACTER*4 CDER(2)
80 : DATA CDER/'DGAU','DDRO'/
81 : LOGICAL LCDER(2)
82 : REAL*8 XCDER(2)
83 : INTEGER OOOVAL
84 :
85 : -INC SMLREEL
86 : POINTEUR MLREE4.MLREEL,MLDERS.MLREEL,MLDER.MLREEL
87 : -INC CCREEL
88 : -INC SMCHPOI
89 : -INC SMELEME
90 : -INC CCOPTIO
91 : -INC SMCHAML
92 : -INC SMTABLE
93 : -INC SMEVOLL
94 : SEGMENT TABR
95 : REAL*8 TEMA(LOR)
96 : ENDSEGMENT
97 :
98 : * syntaxe 3 ?
98 : CALL LIROBJ('NUAGE',INUA,0,IRETOU)
99 : IF(IRETOU.EQ.1) THEN
99 : CALL IPLNU1(INUA)
100 : RETURN
101 : ENDIF
102 : * syntaxe 2 ?
102 : CALL LIROBJ('TABLE',IPOT,0,IRETOU)
103 : IF(IRETOU.EQ.1) THEN
103 : CALL LIRREE(TEMPS,1,IRETOU)
104 : IF(IERR.NE.0) RETURN
105 : GOTO 50
106 : ENDIF
107 :
108 : * syntaxe 1 (INDIC=1 à 4)
108 : CALL LIRREE(TEMPS,0,IRETOU)
109 : IF(IRETOU.EQ.1) THEN
110 : INDIC=1
111 : ELSE
111 : CALL LIROBJ('CHPOINT ',MTEMP,0,IRETOU)
112 : IF(IRETOU.EQ.1) THEN
113 : INDIC=2
114 : ELSE
114 : CALL LIROBJ('MCHAML',IPO1,0,IRETOU)
115 : IF(IRETOU.EQ.1) THEN
116 : INDIC=4
117 : ELSE
117 : CALL LIROBJ('LISTREEL',MLIST,0,IRETOU)
118 : IF(IRETOU.EQ.0) THEN
119 : * Pas d opérande correcte trouvée
119 : CALL QUETYP(MOTERR(1:8),0,IRETOU)
120 : IF(IRETOU.NE.0) THEN
121 : * On ne veut pas d'objet de type %m1:8
121 : CALL ERREUR(39)
122 : ELSE
123 : * Cet opérateur a encore besoin d'un opérande
123 : CALL ERREUR(533)
124 : ENDIF
125 : RETURN
126 : ENDIF
127 : INDIC=3
128 : ENDIF
129 : ENDIF
130 : ENDIF
131 :
132 : * Lecture de la fonction à interpoler
132 : CALL LIROBJ('EVOLUTIO',IEVL,0,IRETOU)
133 : IF(IRETOU.EQ.1) THEN
134 : MEVOLL=IEVL
135 : SEGACT, MEVOLL
136 : * Vérification que l'évolution est élémentaire
137 : IF (IEVOLL(/1).NE.1) THEN
138 : * Opération interdite sur un objet complexe
138 : CALL ERREUR(25)
139 : RETURN
140 : ENDIF
141 : KEVOLL=IEVOLL(1)
142 : SEGACT, KEVOLL
143 : * de type scalaire
144 : IF (ITYEVO.NE.'REEL') THEN
145 : * Opération interdite sur un objet complexe
145 : CALL ERREUR(25)
146 : RETURN
147 : ENDIF
148 : * peuplée de flottants
149 : IF ((TYPX.NE.'LISTREEL').OR.(TYPY.NE.'LISTREEL')) THEN
150 : * Certaines courbes ne sont pas du bon type
150 : CALL ERREUR(871)
151 : RETURN
152 : ENDIF
153 : KTE=IPROGX
154 : KFT=IPROGY
155 : SEGDES,KEVOLL,MEVOLL
156 : ELSE
156 : CALL LIROBJ('LISTREEL',KTE,1,IRETOU)
157 : IF(IERR.NE.0) RETURN
157 : CALL LIROBJ('LISTREEL',KFT,1,IRETOU)
158 : IF(IERR.NE.0) RETURN
159 : ENDIF
160 :
161 : * Longueur des listes décrivant la fonction
162 : MLREE1=KTE
163 : MLREE2=KFT
164 : SEGACT MLREE1,MLREE2
165 : IF(MLREE1.PROG(/1).NE.MLREE2.PROG(/1)) THEN
166 : * Les suites n ont pas la même longueur
166 : CALL ERREUR(212)
167 : RETURN
168 : ENDIF
169 : LON=MLREE1.PROG(/1)
170 : C 897 2
171 : C La dimension des LISTREEL doit etre plus grande que 1
172 : IF (LON.LT.1) THEN
172 : CALL ERREUR(897)
173 : RETURN
174 : ENDIF
175 : TDEB=MLREE1.PROG(1)
176 : TFIN=MLREE1.PROG(LON)
177 : * Les x sont-ils croissants ou décroissants ?
178 : IF (TFIN.GE.TDEB) THEN
179 : ISENS=0
180 : ELSE
181 : ISENS=1
182 : JG=LON
183 : SEGINI,MLREE3
184 : SEGINI,MLREE4
185 : DO ILON=1,LON
186 : MLREE3.PROG(ILON)=MLREE1.PROG(LON-ILON+1)
187 : MLREE4.PROG(ILON)=MLREE2.PROG(LON-ILON+1)
188 : ENDDO
189 : SEGDES,MLREE1
190 : SEGDES,MLREE2
191 : MLREE1=MLREE3
192 : MLREE2=MLREE4
193 : TDEB=MLREE1.PROG(1)
194 : TFIN=MLREE1.PROG(LON)
195 : ENDIF
196 : * On va vérifier que la liste est ordonnée
197 : IVERI=1
198 : C
199 : C Vérification que la liste est ordonnée
200 : C
201 : IF (IVERI.EQ.1) THEN
202 : TPRE=TDEB
203 : DO ILON=2,LON
204 : TCOU=MLREE1.PROG(ILON)
205 : IF (TCOU.LT.TPRE) GOTO 6661
206 : TPRE=TCOU
207 : ENDDO
208 : ENDIF
209 : GOTO 1
210 : 6661 CONTINUE
211 : C "La suite de réels doit être croissante"
212 : C "monotone" !
212 : CALL ERREUR(249)
213 : RETURN
214 : 1 CONTINUE
215 : *
216 : * Y a-t-il une option SPLINE ?
217 : *
217 : CALL LIRMOT(CLE,1,ISPLIN,0)
218 : IF (ISPLIN.EQ.1) THEN
219 : LCDER(1)=.FALSE.
220 : LCDER(2)=.FALSE.
221 : * Lecture des mots clés et valeurs associées
222 : * On lit les valeurs des dérivées premières à gauche et à droite
223 : * Si elles ne sont pas données, c'est la condition à la limite
224 : * naturelle qui s'applique
225 : 77 CONTINUE
225 : CALL LIRMOT(CDER,2,ICDER,0)
226 : IF (ICDER.GT.0) THEN
227 : LCDER(ICDER)=.TRUE.
227 : CALL LIRREE(XCDER(ICDER),1,IRETOU)
228 : IF(IERR.NE.0) RETURN
229 : GOTO 77
230 : ENDIF
231 : JG=LON
232 : SEGINI MLDERS
233 : SEGINI MLDER
234 : IF (LCDER(1)) THEN
235 : * Cas où on prescrit la dérivée première à gauche
236 : MLDERS.PROG(1)=-0.5D0
237 : DX=MLREE1.PROG(2)-MLREE1.PROG(1)
238 : DY=MLREE2.PROG(2)-MLREE2.PROG(1)
239 : MLDER.PROG(1)=(3.D0/DX)*((DY/DX)-XCDER(1))
240 : ELSE
241 : * Condition de bord naturelle (dérivée seconde nulle)
242 : MLDERS.PROG(1)=XZERO
243 : MLDER.PROG(1)=XZERO
244 : ENDIF
245 : DO ILON=2,LON-1
246 : XIM=MLREE1.PROG(ILON-1)
247 : XI=MLREE1.PROG(ILON)
248 : XIP=MLREE1.PROG(ILON+1)
249 : YIM=MLREE2.PROG(ILON-1)
250 : YI=MLREE2.PROG(ILON)
251 : YIP=MLREE2.PROG(ILON+1)
252 : DXIM=XI-XIM
253 : DXI2=XIP-XIM
254 : DXIP=XIP-XI
255 : XRAP=DXIM/DXI2
256 : XP=XRAP*MLDERS.PROG(ILON-1)+2.D0
257 : MLDERS.PROG(ILON)=(XRAP-1.D0)/XP
258 : DYIP=YIP-YI
259 : DYIM=YI-YIM
260 : MLDER.PROG(ILON)=(6.D0*(DYIP/DXIP-DYIM/DXIM)/DXI2-XRAP
261 : $ *MLDER.PROG(ILON-1))/XP
262 : ENDDO
263 : IF (LCDER(2)) THEN
264 : XQN=0.5D0
265 : DX=MLREE1.PROG(LON)-MLREE1.PROG(LON-1)
266 : DY=MLREE2.PROG(LON)-MLREE2.PROG(LON-1)
267 : XUN=(3.D0/DX)*(XCDER(2)-(DY/DX))
268 : ELSE
269 : * Condition de bord naturelle (dérivée seconde nulle)
270 : XQN=0.D0
271 : XUN=0.D0
272 : ENDIF
273 : MLDERS.PROG(LON)=(XUN-XQN*MLDER.PROG(LON-1))/
274 : $ (XQN*MLDERS.PROG(LON-1)+1.D0)
275 : DO ILON=LON-1,1,-1
276 : MLDERS.PROG(ILON)=MLDERS.PROG(ILON)*MLDERS.PROG(ILON+1)
277 : $ +MLDER.PROG(ILON)
278 : ENDDO
279 : SEGSUP MLDER
280 : ELSE
281 : MLDERS=0
282 : ENDIF
283 : * On ne va pas renvoyer d'erreur si t est en-dehors des bornes
284 : * On renvoie la valeur de la fonction pour t le plus proche
285 : IVERI=0
286 : *
287 : * Répartition suivant le type de l'objet fourni
288 : *
289 : GOTO (10,20,30,40) INDIC
290 :
291 : ****************** T0 FLOTTANT *******************************
292 :
293 : 10 CONTINUE
293 : CALL INTER5(TEMPS,MLREE1,MLREE2,FT0,ISPLIN,MLDERS,IVERI,IRET)
293 : IF (IRET.EQ.1) CALL ECRREE(FT0)
294 : GOTO 999
295 :
296 : ********************* T0 CHPOINT *****************************
297 :
298 : 20 CONTINUE
299 : MCHPO1=MTEMP
300 : SEGINI,MCHPOI=MCHPO1
301 : MFT0=MCHPOI
302 : NS=IPCHP(/1)
303 : DO 21 IA=1,NS
304 : MSOUP1=IPCHP(IA)
305 : SEGINI,MSOUPO=MSOUP1
306 : NC=NOHARM(/1)
307 : IPCHP(IA)=MSOUPO
308 : IPT1=IGEOC
309 : SEGINI,IPT2=IPT1
310 : SEGDES IPT2
311 : IGEOC=IPT2
312 : MPOVA1=IPOVAL
313 : SEGINI,MPOVAL=MPOVA1
314 : IPOVAL=MPOVAL
315 : N=VPOCHA(/1)
316 : DO IB=1,N
317 : DO IC=1,NC
318 : TEMPS=VPOCHA(IB,IC)
318 : CALL INTER5(TEMPS,MLREE1,MLREE2,FT0,ISPLIN,MLDERS,
319 : $ IVERI,IRET)
320 : IF (IRET.EQ.0) THEN
321 : SEGSUP MCHPOI,MSOUPO,MPOVAL
322 : GOTO 999
323 : ENDIF
324 : VPOCHA(IB,IC)=FT0
325 : ENDDO
326 : ENDDO
327 : SEGDES MPOVAL
328 : SEGDES MSOUPO
329 : 21 CONTINUE
330 : SEGDES MCHPOI
330 : CALL ECROBJ('CHPOINT ',MFT0)
331 : GOTO 999
332 :
333 : ******************* T0 EST UN LISTREEL ***********************
334 :
335 : 30 CONTINUE
336 : MLREE3=MLIST
337 : SEGACT MLREE3
338 : LONG=MLREE3.PROG(/1)
339 : JG=LONG
340 : SEGINI MLREEL
341 : MSOL=MLREEL
342 : DO 31 ILOOP=1,LONG
343 : TEMPS=MLREE3.PROG(ILOOP)
343 : CALL INTER5(TEMPS,MLREE1,MLREE2,FTO,ISPLIN,MLDERS,
344 : $ IVERI,IRET)
345 : IF(IRET.EQ.0) GOTO 999
346 : PROG(ILOOP)=FTO
347 : 31 CONTINUE
348 : SEGDES MLREEL,MLREE3
348 : CALL ECROBJ('LISTREEL',MSOL)
349 : GOTO 999
350 :
351 : *********************** T0 MCHAML ***************************
352 :
353 : 40 CONTINUE
354 : IRET=0
355 : MCHEL1=IPO1
356 : SEGINI,MCHELM=MCHEL1
357 : IRET=MCHELM
358 : NSOUS=IMACHE(/1)
359 : DO 72 IA=1,NSOUS
360 : MCHAM1=ICHAML(IA)
361 : SEGINI,MCHAML=MCHAM1
362 : ICHAML(IA)=MCHAML
363 : DO 75 ICOMP=1,IELVAL(/1)
364 : MELVA1 = IELVAL(ICOMP)
365 : SEGINI,MELVAL=MELVA1
366 : IELVAL(ICOMP) = MELVAL
367 : SEGACT MELVA1
368 : IF (TYPCHE(ICOMP).NE.'REAL*8') GOTO 166
369 : N1PTEL=VELCHE(/1)
370 : N1EL =VELCHE(/2)
371 : N2PTEL=0
372 : N2EL =0
373 : DO IB=1,N1PTEL
374 : DO ID=1,N1EL
375 : TEMPS=MELVA1.VELCHE(IB,ID)
375 : CALL INTER5(TEMPS,MLREE1,MLREE2,FT0,ISPLIN,
376 : $ MLDERS,IVERI,IREE)
377 : IF (IREE.EQ.0) THEN
378 : SEGSUP MCHELM,MCHAML,MELVAL
379 : SEGDES MELVA1
380 : GOTO 999
381 : ENDIF
382 : VELCHE(IB,ID)=FT0
383 : ENDDO
384 : ENDDO
385 : 166 SEGDES MELVAL,MELVA1
386 : 75 CONTINUE
387 : SEGDES MCHAML
388 : 72 CONTINUE
389 : SEGDES MCHELM
389 : CALL ECROBJ('MCHAML',IRET)
390 : GOTO 999
391 :
392 : ************************ OBJET1 TABLE ******************************
393 :
394 : 50 CONTINUE
395 : MTABLE = IPOT
396 : SEGACT MTABLE
397 : LO = MLOTAB
398 :
399 : *-- Vérification du format de la table
400 : IF (LO.LE.2) THEN
401 : * La table n'a pas le format désiré
401 : CALL ERREUR(647)
402 : SEGDES MTABLE
403 : RETURN
404 : ENDIF
405 :
406 : LOR = LO
407 : SEGINI TABR
408 :
409 : *-- Vérification du sous-type de la table
410 : * IMOT est son indice dans la table
411 : IOK = 0
412 : DO 55 I=1,LO
413 : TYPE1 = MTABTI(I)
414 : IF(TYPE1.EQ.'MOT ') THEN
415 : CHARIN = 'SOUSTYPE'
416 : TYPOBJ = ' '
416 : CALL ACCTAB(MTABLE,TYPE1,IVALIN,XVALIN,CHARIN,LOGIN,
417 : $ IOBIN,TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE)
418 : segact mtable
419 : IF(CHARIN.EQ.'SOUSTYPE') THEN
420 : IF(CHARRE.EQ.'RESULTAT') THEN
421 : IOK = 1
422 : IMOT = I
423 : ENDIF
424 : ENDIF
425 : ENDIF
426 : 55 CONTINUE
427 :
428 : IF(IOK.EQ.0) THEN
429 : * Le sous-type de la table est incorrect
429 : CALL ERREUR(648)
430 : SEGDES MTABLE
431 : SEGSUP TABR
432 : RETURN
433 : ENDIF
434 :
435 : *-- Vérification qu'on a bien des flottants en indice
436 : * en dehors du sous-type
437 : J = 0
438 : DO 56 I=1,LO
439 : IF (IMOT.EQ.I) GOTO 56
440 : J=J+1
441 : TYPEI = MTABTI(I)
442 :
443 : IF (TYPEI.NE.'FLOTTANT') THEN
444 : * La table n'a pas le format desire
444 : CALL ERREUR(647)
445 : SEGDES MTABLE
446 : SEGSUP TABR
447 : RETURN
448 : ENDIF
449 :
450 : TEMA(J) = RMTABI(I)
451 : 56 CONTINUE
452 :
453 : *-- Vérification de l'ordonnancement des indices
454 : DO 57 I=1,LOR-2
455 : TEM1 = TEMA(I)
456 : TEM2 = TEMA(I+1)
457 :
458 : IF(TEM1.GT.TEM2) THEN
459 : * la liste des indices n'est pas ordonnee
459 : CALL ERREUR(211)
460 : SEGDES MTABLE
461 : SEGSUP TABR
462 : RETURN
463 : ENDIF
464 :
465 : 57 CONTINUE
466 :
467 : IF(IMOT.EQ.1) THEN
468 : TEM1 = RMTABI(2)
469 : TYPO1 = MTABTV(2)
470 : IVALO1 = MTABIV(2)
471 : NDEB = 3
472 : ENDIF
473 : IF(IMOT.EQ.2) THEN
474 : TEM1 = RMTABI(1)
475 : TYPO1 = MTABTV(1)
476 : IVALO1 = MTABIV(1)
477 : NDEB = 3
478 : ENDIF
479 : IF(IMOT.GE.3) THEN
480 : TEM1 = RMTABI(1)
481 : TYPO1 = MTABTV(1)
482 : IVALO1 = MTABIV(1)
483 : NDEB = 2
484 : ENDIF
485 :
486 : DO 58 I=NDEB,LOR
487 : IF(IMOT.EQ.I) GOTO 58
488 : TEM2 = RMTABI(I)
489 : TYPO2 = MTABTV(I)
490 : IVALO2 = MTABIV(I)
491 : IF((TEM1.LE.TEMPS).AND.(TEMPS.LE.TEM2)) THEN
492 : DTEM = (TEMPS-TEM1)/(TEM2-TEM1)
493 : IF(TYPO1.EQ.'CHPOINT ') THEN
494 : IF(TYPO2.EQ.'CHPOINT ') THEN
494 : CALL ECROBJ('CHPOINT',IVALO2)
494 : CALL ECROBJ('CHPOINT',IVALO1)
495 : ELSE
495 : CALL ERREUR(647)
496 : SEGDES MTABLE
497 : SEGSUP TABR
498 : RETURN
499 : ENDIF
500 : ENDIF
501 : IF(TYPO1.EQ.'MCHAML ') THEN
502 : IF(TYPO2.EQ.'MCHAML ') THEN
502 : CALL ECROBJ('MCHAML',IVALO2)
502 : CALL ECROBJ('MCHAML',IVALO1)
503 : ELSE
503 : CALL ERREUR(647)
504 : SEGDES MTABLE
505 : SEGSUP TABR
506 : RETURN
507 : ENDIF
508 : ENDIF
508 : CALL ECRREE(DTEM)
508 : CALL ECRREE(1.D0-DTEM)
508 : CALL COLI
509 : SEGDES MTABLE
510 : SEGSUP TABR
511 : GOTO 500
512 : ENDIF
513 : TEM1 = TEM2
514 : TYPO1 = TYPO2
515 : IVALO1 = IVALO2
516 : 58 CONTINUE
517 :
518 : * Le temps est en dehors des limites de la table
518 : CALL ERREUR(210)
519 : SEGDES MTABLE
520 : SEGSUP TABR
521 : GOTO 500
522 :
523 : 999 CONTINUE
524 : IF (ISENS.EQ.0) THEN
525 : SEGDES MLREE1
526 : SEGDES MLREE2
527 : ELSE
528 : SEGSUP MLREE1
529 : SEGSUP MLREE2
530 : ENDIF
531 : IF (ISPLIN.EQ.1) THEN
532 : SEGSUP MLDERS
533 : ENDIF
534 : * Sortie
535 : 500 CONTINUE
536 : RETURN
537 : END
538 :
539 :
540 :
541 :
542 :
543 :
544 :
545 :
546 :
547 :
548 :
© Cast3M 2003 - All rights reserved.
Disclaimer