integral.pro
2.96 KB
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
function integral,tab_1,tab_2,v1,v2,silent=silent,double=double
;+
; NAME:
; integral
; CALLING SEQUENCE:
; res=integral(tab_1,tab_2,v1,v2)
; PURPOSE:
; integrate into tab3(tab2) from v1 to v2
; INPUTS:
; tab_1 = array (1D) of x values
; tab_2 = array (1D) of y values
; v1 = inf value of the integral
; v2 = sup value of the integral
; OPTIONAL INPUT:
; silent = if set routine behaves silently
; OUTPUTS:
; result of the integral
; PROCEDURE AND SUBROUTINE USED
; SIDE EFFECTS:
; None
; MODIFICATION HISTORY:
; written JPB 09-08-92
; modified July 22nd 1994 to handle intervals < point separations
; modified DJM June 10th 2010 to perform calculation in double
;-
IF N_PARAMS(0) EQ 0 THEN begin
print,'Calling Sequence: res=integral(x,y,x1,x2)'
print,'Accepted Key-words: /silent'
goto,sortie
endif
IF keyword_set(double) then half=0.5D else half=0.5
;the aim is to integrate into tab3(tab2) from v1 to v2
;integral in perform linearily
if keyword_set(double) then begin
tab1=double(tab_1)
tab2=double(tab_2)
endif else begin
tab1=tab_1
tab2=tab_2
endelse
si1=size(tab1)
si1=si1(3)
si2=size(tab2)
si2=si2(3)
if si1 ne si2 then begin
print,'pb with array dimensions'
goto,sortie
endif
mi=min(tab1)
ma=max(tab1)
IF keyword_set(double) then coef=1.D else coef=1.
IF tab_1(0) GT tab_1(1) THEN BEGIN
IF n_elements(silent) EQ 0 THEN BEGIN
print,'Array in decreasing order was reordered'
ENDIF
for i=0,si1-1 do tab1(i)=tab_1(si1-i-1)
for i=0,si1-1 do tab2(i)=tab_2(si1-i-1)
ENDIF
if v1 lt mi or v1 gt ma then begin
IF n_elements(silent) EQ 0 THEN BEGIN
message,'first value out of range'
ENDIF
goto,sortie
endif
if v2 lt mi or v2 gt ma then begin
message,'second value out of range'
goto,sortie
endif
if v1 gt v2 then begin
a=v2
v_2=v1
v_1=a
IF keyword_set(double) then coef=-1.D else coef=-1.
endif
if v1 le v2 then begin
v_2=v2
v_1=v1
IF keyword_set(double) then coef=1.D else coef=1.
endif
ind=where(tab1 lt v_2 and tab1 ge v_1,count)
IF count NE 0 THEN BEGIN
a1=tab2(ind)
IF keyword_set(double) then a1(*)=0.D else a1(*)=0.
a1=(tab1(ind+1)-tab1(ind))*(tab2(ind)+half*(tab2(ind+1)-tab2(ind)))
sia1=size(a1)
sia1=sia1(3)
I3=total(a1)-a1(sia1-1)
mi_tab1=min(tab1(ind))
ma_tab1=max(tab1(ind))
mi_tab2=tab2(where(tab1 eq mi_tab1))
ma_tab2=tab2(where(tab1 eq ma_tab1))
IF keyword_set(double) then truc = dblarr(1) else truc=fltarr(1)
truc(0)=v_1
tab2_v1=INTERPOL(TAB2,TAB1,V_1)
truc(0)=v_2
tab2_v2=interpol(tab2,tab1,v_2)
I1=(mi_tab1-v_1)*(tab2_v1+half*(mi_tab2-tab2_v1))
I2=(v_2-ma_tab1)*(ma_tab2+half*(tab2_v2-ma_tab2))
Integ=I1+I2+I3
; stop
ENDIF ELSE BEGIN
ind=where(tab_1 LT v_1) & vv_1=max(tab_1(ind))
; ind=where(tab_1 gT v_2) & vv_2=min(tab_1(ind))
ind=where(tab_1 GE v_2) & vv_2=min(tab_1(ind))
yy_1=tab_2(where(tab_1 EQ vv_1)) & yy_2=tab_2(where(tab_1 EQ vv_2))
a=(yy_2-yy_1)/(vv_2-vv_1) & b=yy_1-a*vv_1
Integ=half*a*(v_2^2-v_1^2)+b*(v_2-v_1)
ENDELSE
return,integ
sortie:
end