Discussion :

[noomane]

Bonjour tout le monde, j’ai le problème suivant : je veux faire l’intégrale suivante :


Avec est la fonction suivante :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
!LA FONCTION DENSITE DE PROBABILITE DE LA LOI NORMALE  
!----------------------------------------------------
 
	FUNCTION DNORM(T)
	INCLUDE 'comm.f'
 
	PI=3.141592654
 
	DNORM=(1/(EC*SQRT(2*PI)))*EXP(-0.5*(((T-AVG)/EC)**2))
 
	RETURN
 
	END

EC et AVG sont des constantes avec EC=900.0 et AVG=90.0

Et la fonction (.) est la fonction avec f est même que avec les paramètre EC1=80.0 et AVG1=850.0

Donc on aura besoin de la fonction suivante :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
!LA FONCTION DENSITE DE PROBABILITE DE LA LOI NORMALE  
!----------------------------------------------------
 
	FUNCTION DNORM1(T)
	INCLUDE 'comm.f'
 
	PI=3.141592654
 
	DNORM1=(1/(EC1*SQRT(2*PI)))*EXP(-0.5*(((T-AVG1)/EC1)**2))
 
	RETURN
 
	END

Pour faire l’intégrale on besoin d’une subroutine d’intégrale QTRAP :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
!CALCUL D'UNE INTEGRALE PAR LA METHODE DU TRAPEZE
!------------------------------------------------
 
	SUBROUTINE QTRAP(FUNC,U,V,S)
	PARAMETER (EPS=1.E-2,JMAX=20)
	EXTERNAL FUNC
	IF ((V-U).LT.1.E-4) THEN
 
	S=1.E-6
 
	RETURN
	ENDIF
 
	IF (V.EQ.0.0) THEN
 
	V=0.001
 
	ENDIF
 
	OLDS=-1.E30
 
	DO 10 J=1,JMAX
 
	 CALL TRAPZD(FUNC,U,V,S,J)
 
	 IF (ABS(S-OLDS).LT.EPS*ABS(OLDS))RETURN
 
	 OLDS=S
 
 10	CONTINUE
 
	 END
 
!SUITE DE LA METHODE DU TRAPEZE
!------------------------------
 
	SUBROUTINE TRAPZD(FUNC,U,V,S,N)
	EXTERNAL FUNC
	IF (N.EQ.1)THEN
	 S=0.5*(V-U)*(FUNC(U)+FUNC(V))
	IT=1
 
	ELSE
	 TNM=IT
	 DEL=(V-U)/TNM
	 G=U+0.5*DEL
	 SUM=0.
 
	 DO 12 J=1,IT
	  SUM=SUM+FUNC(G)
	  G=G+DEL
 
 12	  CONTINUE
	 S=0.5*(S+(V-U)*SUM/TNM)
	 IT=2*IT
 
	ENDIF
	RETURN
	END

Donc la fonction CNORM qui est F = est représentée par la fonction suivante :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
! LA FONCTION DE DISTRIBUTION DE LA LOI NORMALE
!----------------------------------------------
 
	FUNCTION CNORM(T)
	INCLUDE 'comm.f'
	EXTERNAL QTRAP,DNORM1
 
	CALL QTRAP(DNORM,0.0,T,S)
 
	CNORM=S
 
	RETURN
 
	END

Le problème c’est toujours la même valeur que je trouve qui’ est 1.0000E-0.6

Comment faire ?

Merci.

[genteur slayer]

où se trouve l'erreur à la sortie de CNORM ou de la subroutine c'est quoi qui rentre dedans c quoi qui sort....
l'inplicit c'est partique mais ça cause des erreurs, de plus utilise les balises de code, ton msg sera plus clair. (autre petit conseil en passant: lache le fortran77, passe à un plus récent 90/95 ou 2003)

au fait c quoi l'intégrale que tu veux calculer? tu veut l'espérence de ta loi normale??? si c le cas, t'embête pas colle direct le résultat...

[noomane]

merci pour ton aide, j'ai fait joindre un fichier .doc ou j'explique tout, en attente d'une reponse rapide,

noomane.

[genteur slayer]

petites astuces au passage, pour la constante Pi, plutôt que de la rentrer soit même avec des précision aléeatoire tu peut donner Pi=4.d0*datan(1.d0) comme ça la machine te donne la précision maximal dont elle est capable.
en deux: commente à fond!!!! avant chaque fonctions tu écris en français ce qu'elle fait (du coup on voit très vite quand il y a un désaccord avec le language) en fortran, le "!" est équivalent au "//" du c/c++ tu peut le mettre n'importe où, la fin de la ligne est commentaire... si tu utilise le format fixe de fortran, tu peut écrire après la 72ième colonne (ou la 132 selon comment tu regle ton compilo) c'est des commentaire... (même si papa GNU dit qu'un beaux prog n'a pas besoin de commentaire pour être compréhenssible, je te conseil ce ratio: une ligne de code avec 2 lignes de commentaire (en moyenne))

pour trouvé ton erreur, je te conseil de tester toute les fonctions et subroutine une par une. puis de regarder la visibilité de tes variables (t'aurais pas une variable locale à qui tu donne le resultats et se perd dans la mémoire) je suis pas du tout rompu à l'utilisation du common (à l'école on m'a tjrs dit de les éviter) mais ça peut venir de là... tu devrais aussi verifié l'orthographe de tes variables (comme tu es en implicit pas d'erreur de compil ni d'exec) condeil donc: calle un "implicit none" au début de chaque fonction et déclare explicitement toute tes variable (je sais c long et chiant mais le débug l'est encore plus)

sinon je ne connais pas du tout ta version de fortran, elle tourne sur quoi?

[noomane]

Merci à vous, j'ai résolu le prob mais ça prend beaucoup de temps de calcul malgré que ma machine soit performante.

Voici le code :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
		! prog conditionnelle
include 'comm.f'
external chvar,DNORM,CNORM,QTRAP1,TRAPZD1,QTRAP2,TRAPZD2,QTRAP3,TRAPZD3,QTRAP4,TRAPZD4,DNORM1,CMTN,zbrent,CCTM,CMTN1
real XI(100,100)
real PS(100)
!real KS
AVG=900.0
EC=90.0
AVG1=850.0
EC1=80.0
H=10.0
PS(1)=0.1
B=785.0
j=1.0
 do while (PS(j).LT.1.0)
 
PK=0.0
 
XI(0,j)=0.0
	do i=1,2
 
		!determiantion PHI(Xi-1)
 
		call QTRAP3(DNORM1,0.0,XI(i-1,j),ST)
 
		!calcul de la quantité PK=Ps+(Ps-1)+PHI(Xi-1):
 
		PK=PS(j)-(PS(j)-1)*ST
 
		!annulation de la fonction CCTM pour determiner les Xi:
 
		XI(i,j)=zbrent(CCTM,-10000.0,10000.0,0.0001)
 
 
 
 
		!C=XI(i,j)
 
		!write(6,*)'valeur XI+1',C
 
 
		write(6,*)'X',i,XI(i,j)
 
		!call QTRAP2(CMTN,0.0,B,S)
 
		!pk2=S
 
		!write(6,*)'val integ',pk2
 
	end do
			do i=0,2
 
				B=XI(i,j)
 
				call QTRAP2(CMTN,0.0,B,S)
 
				pk2=S
 
				write(6,*)'val integ',pk2
 
				C=XI(i+1,j)
 
				call QTRAP4(CMTN1,0.0,C,S1)
 
				pk3=S1
 
				write(6,*)'val integ1',pk3
 
			end do
 
write(6,*)'la sequence optimal avec P=',PS(j)
 
j=j+1
 
PS(j)=PS(j-1)+0.1
 
 end do
 
 
end
 
! la fonction qui determine les xi : CCTM = PHI(Xi)-Ps+(Ps-1)+PHI(Xi-1) 
 
	FUNCTION CCTM(T)
	include 'comm.f'
	external CNORM
 
	CCTM = CNORM(T) - PK
 
	return
	end
 
 
 ! LA DERIVE DU COUT MOYEN TOTAL PAR UNITE DE TEMPS POUR LA LOI NORMALE
!---------------------------------------------------------------------
 
	FUNCTION CMTN1(T)
	INCLUDE 'comm.f'
	EXTERNAL CNORM,DNORM1
 
	CMTN1= DNORM1(T) * CNORM(C-T)
	!write(6,*)'valeur de T1',t
 
	RETURN
 
	END
 
 
! LA DERIVE DU COUT MOYEN TOTAL PAR UNITE DE TEMPS POUR LA LOI NORMALE
!---------------------------------------------------------------------
 
	FUNCTION CMTN(T)
	INCLUDE 'comm.f'
	EXTERNAL CNORM,DNORM1
 
	CMTN= DNORM1(T) * CNORM(B-T)
	write(6,*)'valeur de T',t
 
	RETURN
 
	END
 
 !LA FONCTION DENSITE DE PROBABILITE DE LA LOI NORMALE 
!----------------------------------------------------
 
	FUNCTION DNORM1(T)
	INCLUDE 'comm.f'
 
	PI=3.141592654
 
	DNORM1=(1/(EC1*SQRT(2*PI)))*EXP(-0.5*(((T-AVG1)/EC1)**2))
 
	RETURN
 
	END
 
 
 !LA FONCTION DENSITE DE PROBABILITE DE LA LOI NORMALE 
!----------------------------------------------------
 
	FUNCTION DNORM(T)
	INCLUDE 'comm.f'
 
	PI=3.141592654
 
	DNORM=(1/(EC*SQRT(2*PI)))*EXP(-0.5*(((T-AVG)/EC)**2))
 
	RETURN
 
	END
 
  ! LA FONCTION DE DISTRIBUTION DE LA LOI NORMALE
!----------------------------------------------
 
	FUNCTION CNORM(T)
	INCLUDE 'comm.f'
	EXTERNAL QTRAP1,DNORM
 
	CALL QTRAP1(DNORM,0.0,T,S)
 
	CNORM=S
 
	!write(6,*)'val de F(XI-T)',S
 
	RETURN
 
	END
 
 !CALCUL D'UNE INTEGRALE PAR LA METHODE DU TRAPEZE
!------------------------------------------------
 
	SUBROUTINE QTRAP1(FUNC,U,V,S)
	PARAMETER (EPS=1.E-2,JMAX=20)
	EXTERNAL FUNC
	IF ((V-U).LT.1.E-4) THEN
 
	S=1.E-6
 
	RETURN
	ENDIF
 
	IF (V.EQ.0.0) THEN
 
	V=0.001
 
	ENDIF
 
	OLDS=-1.E30
 
	DO 10 J=1,JMAX
 
	 CALL TRAPZD1(FUNC,U,V,S,J)
 
	 IF (ABS(S-OLDS).LT.EPS*ABS(OLDS))RETURN
 
	 OLDS=S
 
 10	CONTINUE
 
	 END
 
 !SUITE DE LA METHODE DU TRAPEZE
!------------------------------
 
	SUBROUTINE TRAPZD1(FUNC,U,V,S,N)
	EXTERNAL FUNC
	IF (N.EQ.1)THEN
	 S=0.5*(V-U)*(FUNC(U)+FUNC(V))
	IT=1
 
	ELSE
	 TNM=IT
	 DEL=(V-U)/TNM
	 G=U+0.5*DEL
	 SUM=0.
 
	 DO 12 J=1,IT
	  SUM=SUM+FUNC(G)
	  G=G+DEL
 
 12	  CONTINUE
	 S=0.5*(S+(V-U)*SUM/TNM)
	 IT=2*IT
 
	ENDIF
	RETURN
	END
 
! routine chvar
subroutine chvar(u,v,s)
real u,v
s=u-v
continue 
end
 
!CALCUL D'UNE INTEGRALE PAR LA METHODE DU TRAPEZE
!------------------------------------------------
 
	SUBROUTINE QTRAP2(FUNC,U,V,S)
	PARAMETER (EPS=1.E-2,JMAX=20)
	EXTERNAL FUNC
	IF ((V-U).LT.1.E-4) THEN
 
	S=1.E-6
 
	RETURN
	ENDIF
 
	IF (V.EQ.0.0) THEN
 
	V=0.001
 
	ENDIF
 
	OLDS=-1.E30
 
	DO 10 J=1,JMAX
 
	 CALL TRAPZD2(FUNC,U,V,S,J)
 
	 IF (ABS(S-OLDS).LT.EPS*ABS(OLDS))RETURN
 
	 OLDS=S
 
 10	CONTINUE
 
	 END
 
 !SUITE DE LA METHODE DU TRAPEZE
!------------------------------
 
	SUBROUTINE TRAPZD2(FUNC,U,V,S,N)
	EXTERNAL FUNC
	IF (N.EQ.1)THEN
	 S=0.5*(V-U)*(FUNC(U)+FUNC(V))
	IT=1
 
	ELSE
	 TNM=IT
	 DEL=(V-U)/TNM
	 G=U+0.5*DEL
	 SUM=0.
 
	 DO 12 J=1,IT
	  SUM=SUM+FUNC(G)
	  G=G+DEL
 
 12	  CONTINUE
	 S=0.5*(S+(V-U)*SUM/TNM)
	 IT=2*IT
 
	ENDIF
	RETURN
	END
 
!CALCUL D'UNE INTEGRALE PAR LA METHODE DU TRAPEZE
!------------------------------------------------
 
	SUBROUTINE QTRAP3(FUNC,U,V,S)
	PARAMETER (EPS=1.E-2,JMAX=20)
	EXTERNAL FUNC
	IF ((V-U).LT.1.E-4) THEN
 
	S=1.E-6
 
	RETURN
	ENDIF
 
	IF (V.EQ.0.0) THEN
 
	V=0.001
 
	ENDIF
 
	OLDS=-1.E30
 
	DO 10 J=1,JMAX
 
	 CALL TRAPZD3(FUNC,U,V,S,J)
 
	 IF (ABS(S-OLDS).LT.EPS*ABS(OLDS))RETURN
 
	 OLDS=S
 
 10	CONTINUE
 
	 END
 
 !SUITE DE LA METHODE DU TRAPEZE
!------------------------------
 
	SUBROUTINE TRAPZD3(FUNC,U,V,S,N)
	EXTERNAL FUNC
	IF (N.EQ.1)THEN
	 S=0.5*(V-U)*(FUNC(U)+FUNC(V))
	IT=1
 
	ELSE
	 TNM=IT
	 DEL=(V-U)/TNM
	 G=U+0.5*DEL
	 SUM=0.
 
	 DO 12 J=1,IT
	  SUM=SUM+FUNC(G)
	  G=G+DEL
 
 12	  CONTINUE
	 S=0.5*(S+(V-U)*SUM/TNM)
	 IT=2*IT
 
	ENDIF
	RETURN
	END
 
 
!CALCUL D'UNE INTEGRALE PAR LA METHODE DU TRAPEZE
!------------------------------------------------
 
	SUBROUTINE QTRAP4(FUNC,U,V,S)
	PARAMETER (EPS=1.E-2,JMAX=20)
	EXTERNAL FUNC
	IF ((V-U).LT.1.E-4) THEN
 
	S=1.E-6
 
	RETURN
	ENDIF
 
	IF (V.EQ.0.0) THEN
 
	V=0.001
 
	ENDIF
 
	OLDS=-1.E30
 
	DO 10 J=1,JMAX
 
	 CALL TRAPZD4(FUNC,U,V,S,J)
 
	 IF (ABS(S-OLDS).LT.EPS*ABS(OLDS))RETURN
 
	 OLDS=S
 
 10	CONTINUE
 
	 END
 
 !SUITE DE LA METHODE DU TRAPEZE
!------------------------------
 
	SUBROUTINE TRAPZD4(FUNC,U,V,S,N)
	EXTERNAL FUNC
	IF (N.EQ.1)THEN
	 S=0.5*(V-U)*(FUNC(U)+FUNC(V))
	IT=1
 
	ELSE
	 TNM=IT
	 DEL=(V-U)/TNM
	 G=U+0.5*DEL
	 SUM=0.
 
	 DO 12 J=1,IT
	  SUM=SUM+FUNC(G)
	  G=G+DEL
 
 12	  CONTINUE
	 S=0.5*(S+(V-U)*SUM/TNM)
	 IT=2*IT
 
	ENDIF
	RETURN
	END
 
 
 
 
 
 
 
 
!ALGORITHME DE RESOLUTION D'UNE EQUATION NON LINEAIRE ZBRENT:f(x)=0
!------------------------------------------------------------------
 
	FUNCTION zbrent(func,x1,x2,tol)
    EXTERNAL func
    PARAMETER (ITMAX=100,EPS=3.e-8)
 
    a=x1
    b=x2
    fa=func(a)
	fb=func(b)
	if(fb*fa.GT.0.)pause'Root must be bracketed for zbrent.'
 
    fc=fb
    do 14 iter=1,ITMAX
      if(fb*fc.GT.0.)then
       c=a
       fc=fa
	   d=b-a
       e=d
      endif
	  if(abs(fc).lt.abs(fb)) then
        a=b
        b=c
        c=a
        fa=fb
        fb=fc
        fc=fa
      endif
      tol1=2.*EPS*abs(b)+0.5*tol
      xm=.5*(c-b)
      if(abs(xm).le.tol1 .or. fb.eq.0.)then
       zbrent=b
       return
      endif
      if(abs(e).ge.tol1 .and. abs(fa).gt.abs(fb)) then
       s=fb/fa
       if(a.eq.c) then
         p=2.*xm*s
         q=1.-s
       else
	     q=fa/fc
         r=fb/fc
         p=s*(2.*xm*q*(q-r)-(b-a)*(r-1.))
         q=(q-1.)*(r-1.)*(s-1.)
       endif
       if(p.gt.0.) q=-q
        p=abs(p)
       if(2.*p .lt. min(3.*xm*q-abs(tol1*q),abs(e*q))) then
        e=d
        d=p/q
       else
        d=xm
        e=d
       endif
      else
       d=xm
       e=d
      endif
      a=b
      fa=fb
      if(abs(d) .gt. tol1) then
       b=b+d
      else
	   b=b+sign(tol1,xm)
      endif
      fb=func(b)
14    continue
      pause 'zbrent exceeding maximum iterations'
      zbrent=b
      return
 
	  END

et voici le fichier comm.f :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

	common/param/u,v,s,EC,AVG,EC1,AVG1,B,H,PK,KS,C,S1

essaye de le faire marcher sur votre station, si tu a une autre methode pour le faire tu peux me l'envoyer, et encore merci

[noomane]

g oublier un truc et que si tu peut menvoyer un exemple de fonction ou on utilise l'option implicite none et on declare les variable dans la fonction directement, si tu peut le faire pour mon exemple ca sera tres sympas de votre part, encore une faveur si il existe une mthode de determiner la primitive d'une fonction sur fortran, et merci.

noomane

[genteur slayer]

y a plusieur styles de déclaration, je te donne celle que moi j'utilise (c à dire du fortran 90/95 (je sais pas bien lequel des deux... tantpis)) et je l'applique à une fonction Ray (celle ci définit un profil de rayon pour un tube.... )

le paramètre ki est un entier qui dit combien d'octet prend mon nombre réel, placer ki=kind(1.d0) permet d'avoir des doubles precision sur n'importe quel machine (c le compilo qui s'adapte...) on donne le type du résultat à la manière du C c'est à dire avant d'écrire "function" ensuite l'attirbut "intent' peut avoir 3 valeur :in, out, inout si je détail pas, tu vois à quoi ça sert

ensuite dans la seconde fonction, on voit un bloc interface... là c'est un peut long à expliquer, faut que tu aille voir un cours de fortran 95 disont que sans cela (et avec l'inplicit none) ma seconde fonction ne sait pas du tout ce qu'est 'dsigne' donc on lui dit et on ecopie l'en-tête de fonction (c'est à dire toute les variables présentes sur la première ligne (ici, le x en fait partie donc on doit dire ce que c'est... (j'utilise pas 'ki' à cet endrois car j'avais des pbs de visibilité de la variable et pas envie de me cassé la tête....)

cependant, moi je n'utilise jamais les common, data, et autres donc y a plein d'autres méthodes

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
!=======================================================================
! fonction dsigne(x)
!=======================================================================
real(kind=kind(1.d0)) function dsigne(x)
  implicit none
  integer, parameter :: ki=kind(1.d0)                                  !dimenssionnement des rééls
  real(ki), intent(in) :: x
  dsigne=0.d0
  if (x>0) dsigne=1.d0
  if (x<0) dsigne=-1.d0
  return
end function dsigne
!=======================================================================
 
!=======================================================================
! fonction de rayon initial
!=======================================================================
real(kind=kind(1.d0)) function Ray(x,R,a,omega,mode)
  implicit none
  integer, parameter :: ki=kind(1.d0)                                  !dimenssionnement des rééls
! mode: type de courbe: type 0 constant égal à R
!       type 1 belle sinusoïde moyenne R, amplitude a,
!              nombre de pic par 1 en x omega
!       type 2 vague à base R, crète a, nombre de pic par 1 en x omega
!       type 3 dent de scie (à faire)
!       type 4 créneau même param que sinusoïde
  real(ki) ::  Pi
  real(ki), intent(in)  :: x,a,omega,R
  integer,intent(in):: mode
!-----------------------------------------------------------------------
! bloc interface________________________________________________________
interface
  real(kind=kind(1.d0)) function dsigne(x)
  real(kind=kind(1.d0)), intent(in) :: x
end function dsigne
end interface
!_______________________________________________________________________
  Pi=4.d0*datan(1.d0)
  if (mode==1)then														! une simple sinusoïde paramétrée
	Ray=R+a*dsin(2.d0*Pi*omega*x)
  else if (mode==2)then													! une jolie petite vague
	ray=R-a/2.d0*(0.75*dsin(omega*Pi*x)**6+1.25*dexp(-((dsin(2.d0*Pi*(omega*x-0.125))-1.d0)**2)/0.25))
  else if (mode==4)then													! un crénau
	ray=R+a*dsigne(dsin(2.d0*Pi*omega*x))
  else																	! rayon constant
	Ray=R
  endif
  return
end function ray
!=======================================================================