Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

section.pro @ 76

Last change on this file since 76 was 29, checked in by pinsard, 18 years ago
upgrade of TRIANGULATION according to cerbere.lodyc.jussieu.fr: /usr/home/smasson/SAXO_RD/ : files
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 13.7 KB

Line
1	;------------------------------------------------------------
2	;------------------------------------------------------------
3	;------------------------------------------------------------
4	;+
5	; NAME:
6	;
7	; PURPOSE:
8	;
9	; CATEGORY:
10	;
11	; CALLING SEQUENCE:
12	;
13	; INPUTS:
14	;
15	; KEYWORD PARAMETERS:
16	;
17	; OUTPUTS:
18	;
19	; COMMON BLOCKS:common.pro
20	;
21	; SIDE EFFECTS:
22	;
23	; RESTRICTIONS:
24	;
25	; EXAMPLE:
26	;
27	; MODIFICATION HISTORY:Sebastien Masson (smasson@lodyc.jussieu.fr)
28	;
29	;-
30	;------------------------------------------------------------
31	;------------------------------------------------------------
32	;------------------------------------------------------------
33
34	PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints $
35	, BOXZOOM = boxzoom, TYPE = type, WDEPTH = wdepth $
36	, DIREC = direc, SHOWBUILD = showbuild, ONLYBOX = onlybox $
37	, _extra = ex
38	;
39	;---------------------------------------------------------
40	; include common
41	@cm_4mesh
42	@cm_4data
43	@cm_4cal
44	IF NOT keyword_set(key_forgetold) THEN BEGIN
45	@updatenew
46	@updatekwd
47	ENDIF
48	;--------------------------------------------------------------
49	;---------------------
50	;---------------------
51	; definition de boxzoom en fonction de endpoints
52	; puis redefinition du domaine
53	boxzoom2d = [min([endpoints[0], endpoints[2]], max = ma02), ma02 $
54	, min([endpoints[1], endpoints[3]], max = ma13), ma13]
55	;
56	minprof = 0.
57	profdefault = 200.
58	if n_elements(type) EQ 0 then type = 'nothing'
59	Case N_Elements(Boxzoom) OF
60	0:localbox = [boxzoom2d, minprof, profdefault]
61	1:localbox = [boxzoom2d, minprof, boxzoom[0]]
62	2:localbox = [boxzoom2d, boxzoom[0]]
63	4:if strpos(type, 'z') NE -1 THEN $
64	localbox = [boxzoom2d, minprof, profdefault] ELSE localbox = boxzoom2d
65	5:localbox = [boxzoom2d, minprof, boxzoom[4]]
66	6:localbox = [boxzoom2d, boxzoom[4:5]]
67	Else:BEGIN
68	print, report('Bad definition of the box')
69	stop
70	END
71	ENDCASE
72	nelbox = n_elements(localbox)
73	;
74	if keyword_set(wdepth) then grillechoice = [vargrid, 'W'] $
75	ELSE grillechoice = vargrid
76	domdef, localbox, GRIDTYPE = grillechoice, /findalways, _extra = ex
77	grille, -1, -1, -1, -1, nx, ny
78	; if less than 10 points where found, we apply domdef over the whole domain
79	; -> problem... why 10 points as a test value???
80	; how can we find a good test value?
81	IF nx * ny LE 10 THEN domdef, GRIDTYPE = grillechoice, _extra = ex
82	; on redefinit lon1, ... au cas ou findalways ait ete utilise ds domdef
83	lon1 = min([endpoints[0], endpoints[2]], max = lon2)
84	lat1 = min([endpoints[1], endpoints[3]], max = lat2)
85	; we extend the box along the z axis -> i that way the plot will be drawn
86	; until its bottom part.
87	if strpos(type, 'z') NE -1 THEN BEGIN
88	;on garde les yranges (axe z) avant de changer la boxzoom.
89	!y.range = [localbox[nelbox-1], localbox[nelbox-2]]
90	if vargrid EQ 'W' OR keyword_set(wdepth) then BEGIN
91	firstzw = 0 > (firstzw-1)
92	lastzw = (lastzw+1) < (jpk-1)
93	nzw = lastzw - firstzw + 1
94	ENDIF ELSE BEGIN
95	firstzt = 0 > (firstzt-1)
96	lastzt = (lastzt+1) < (jpk-1)
97	nzt = lastzt - firstzt + 1
98	ENDELSE
99	IF NOT keyword_set(key_forgetold) THEN BEGIN
100	@updateold
101	ENDIF
102	ENDIF
103	!y.range = [localbox[nelbox-1], localbox[nelbox-2]]
104	; on recupere la grille sur la boxzoom
105	;
106	IF NOT keyword_set(onlybox) THEN saveboxparam, 'boxparam4section.dat'
107	grille, -1, -1, -1, -1, nx, ny, nz $
108	, firstx, firsty, firstz, lastx, lasty, lastz
109	; the extend the box in the x and y direction in order to maximise
110	; the number of triangles used to build the section
111	firstx = 0 > (firstx - 1)
112	lastx = (lastx + 1) < (jpi -1)
113	firsty = 0 > (firsty - 1)
114	lasty = (lasty + 1) < (jpj -1)
115
116	domdef, firstx, lastx, firsty, lasty, firstz, lastz $
117	, /index, gridtype = vargrid
118
119	IF keyword_set(onlybox) THEN return
120	;
121	grille, mask, glam, gphi, gdep, nx, ny, nz $
122	, firstx, firsty, firstz, lastx, lasty, lastz
123
124	;---------------------
125	; on definit la triangulation qui va nous permetre de determiner la
126	; section. on la recalcule car elle doit etre definie sur la terre
127	; aussi bien que sur la mer. suivant le sens de la section -plutot
128	; longitude ou plutot latitude- on definit la facon de trianguler
129	if strpos(type, 'x') NE -1 then BEGIN
130	downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2]
131	tri = definetri(nx, ny, (downward)[*])
132	ENDIF ELSE tri = definetri(nx, ny)
133	; If we have an irregular grid that is periodic, then it is possible that
134	; some of the triangle have a very large size (neighborg points on the
135	; sphere but far away when doing the projection) and should not be
136	; taken into account..
137	IF keyword_set(key_irregular) AND keyword_set(key_periodic) THEN BEGIN
138	glamtri = glam[tri]
139	glamtri = abs(glamtri - shift(glamtri, 1, 0))
140	good = temporary(glamtri) LT (10.*max(glam)/nx)
141	good = where(total(temporary(good), 1) EQ 3)
142	tri = (temporary(tri))[*, temporary(good)]
143	ENDIF
144	;---------------------
145	; equation de la droite suivant laquelle on fait la section
146	;---------------------
147	abc = linearequation(endpoints[0:1], endpoints[2:3])
148	glamtri = glam[tri]
149	gphitri = gphi[tri]
150	; quels sont les points de la triangulation qui sont au dessus et au
151	; dessous de la droite ?
152	if abc[1] NE 0 THEN $
153	test = temporary(gphitri) GE (-abc[0]/abc[1]*temporary(glamtri)-abc[2]/abc[1]) $
154	ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0])
155
156	zero123 = total(test, 1)
157	; to keep: triangles de la triangulation qui sont a cheval sur la droite.
158	tokeep1 = where(zero123 EQ 1)
159	tokeep2 = where(temporary(zero123) EQ 2)
160	tokeep = [tokeep1, tokeep2]
161
162	test = test[*, tokeep]
163	tri = tri[*, tokeep]
164	; quel est le sommet du triangle qui est seul d''un cote de la droite?
165	single1 = where(test[*, 0:n_elements(tokeep1)-1] EQ 1)
166	single1 = single1-(single1/3)*3
167	single2 = where(test[*, n_elements(tokeep1):n_elements(tokeep)-1] EQ 0)
168	single2 = single2-(single2/3)*3
169
170	undefine, tokeep
171	undefine, tokeep1
172	undefine, tokeep2
173	undefine, test
174
175	single = [temporary(single1), temporary(single2)]
176	; points1 le point du triangles qui est seul d''un cote de la droite.
177	; point2 l''autre point du triangle de l''autre cote de la droite
178	point1 = [single, single]
179	point2 = [single EQ 0, 1 + (single LE 1)]
180
181	undefine, single
182
183	ntri = (size(tri))[2]
184	index = [lindgen(ntri), lindgen(ntri)]
185
186	points1 = tri[point1, index]
187	points2 = tri[point2, temporary(index)]
188	; points : complexe contenant les couples de points de part et
189	; d''autre de la droite. Ils faut supprimer les doublons.
190	points = dcomplex(points1, points2)
191	points = points[uniq(points, sort(points))]
192	symetrique = dcomplex(imaginary(points), double(points))
193	points = points[where(points-shift(temporary(symetrique), 1) NE 0)]
194	; points1 les coordnnees du point du triangles qui est seul d''un cote de la droite.
195	; point2 les coordnnees de l''autre point du triangle de l''autre cote de la droite
196	points1 = complex(glam[ double(points)], gphi[ double(points)])
197	points2 = complex(glam[imaginary(points)], gphi[imaginary(points)])
198	; droites les equations des droites dont on cherche l''intersection
199	; avec la section.
200	droites = linearequation(points1, points2)
201	inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) )
202
203	; les ccordonnes geographiques des points que l''on cherche sur la section.
204	glamaxe = float(inter)
205	gphiaxe = imaginary(inter)
206	; on les range ds l''ordre croissant entre les bornes de la section
207	if strpos(type, 'x') NE -1 then BEGIN
208	sort = sort(glamaxe)
209	glamaxe = glamaxe[sort]
210	inbox = where(glamaxe GE lon1 AND glamaxe LE lon2)
211	glamaxe = glamaxe[inbox]
212	sort = sort[inbox]
213	gphiaxe = gphiaxe[sort]
214	ENDIF ELSE BEGIN
215	sort = sort(gphiaxe)
216	gphiaxe = gphiaxe[sort]
217	inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2)
218	gphiaxe = gphiaxe[inbox]
219	sort = sort[inbox]
220	glamaxe = glamaxe[sort]
221	ENDELSE
222	points = points[sort]
223	points1 = points1[sort]
224	points2 = points2[sort]
225	inter = inter[sort]
226	poids = abs(points2-inter)/abs(points2-points1)
227
228	;---------------------
229	array = litchamp(field)
230	array = fitintobox(array)
231	if array[0] EQ -1 THEN BEGIN
232	res = -1
233	return
234	ENDIF
235	if n_elements(valmask) EQ 0 THEN valmask = 1e20
236	taille = size(array)
237	if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN
238	direc = 't'
239	array = grossemoyenne(array, 't')
240	taille = size(array)
241	jpt = 1
242	ENDIF
243	case 1 of
244	;----------------------------------------------------------------------------
245	;xy
246	;----------------------------------------------------------------------------
247	taille[0] EQ 2:BEGIN
248	value1 = array[double(points)]
249	terre = where(value1 GT valmask/10)
250	if terre[0] NE -1 then value1[terre] = !values.f_nan
251	value2 = array[imaginary(points)]
252	terre = where(value2 GT valmask/10)
253	if terre[0] NE -1 then value2[terre] = !values.f_nan
254	res = poidsvalue1+(1-poids)value2
255	END
256	;----------------------------------------------------------------------------
257	;xyz
258	;----------------------------------------------------------------------------
259	taille[0] EQ 3 AND jpt EQ 1:BEGIN
260	npoints = n_elements(points)
261	index = double(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz)
262	value1 = array[index]
263	terre = where(value1 GT valmask/10)
264	if terre[0] NE -1 then value1[terre] = !values.f_nan
265	index = imaginary(points)#replicate(1, nz)+replicate(nx*ny, npoints)#lindgen(nz)
266	value2 = array[index]
267	terre = where(value2 GT valmask/10)
268	if terre[0] NE -1 then value2[terre] = !values.f_nan
269	poids = poids#replicate(1, nz)
270	res = poidsvalue1+(1-poids)value2
271	; moyenne suivant z ?
272	if strpos(type, 'z') EQ -1 then begin
273	nan = where(finite(res) EQ 0)
274	if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt]
275	weight = replicate(1, npoints)#e3
276	if nan[0] NE -1 then weight[nan] = !values.f_nan
277	totalweight = total(weight, 2, /nan)
278	zero = where(totalweight EQ 0)
279	if zero[0] NE -1 then totalweight[zero] = !values.f_nan
280	res = total(res*weight, 2, /nan)/totalweight
281	direc = 'z'+string(byte(testvar(var = toto)))
282	endif
283	END
284	;----------------------------------------------------------------------------
285	;xyt
286	;----------------------------------------------------------------------------
287	taille[0] EQ 3 AND jpt NE 1:BEGIN
288	npoints = n_elements(points)
289	index = double(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt)
290	value1 = array[index]
291	terre = where(value1 GT valmask/10)
292	if terre[0] NE -1 then value1[terre] = !values.f_nan
293	index = imaginary(points)#replicate(1, jpt)+replicate(nx*ny, npoints)#lindgen(jpt)
294	value2 = array[index]
295	terre = where(value2 GT valmask/10)
296	if terre[0] NE -1 then value2[terre] = !values.f_nan
297	poids = poids#replicate(1, jpt)
298	res = poidsvalue1+(1-poids)value2
299	END
300	;----------------------------------------------------------------------------
301	;xyzt
302	;----------------------------------------------------------------------------
303	taille[0] EQ 4:BEGIN
304	npoints = n_elements(points)
305	index = double(points)#replicate(1, nzjpt)+replicate(nxny, npoints)#lindgen(nz*jpt)
306	index = reform(index, npoints, nz, jpt, /over)
307	value1 = array[index]
308	terre = where(value1 GT valmask/10)
309	if terre[0] NE -1 then value1[terre] = !values.f_nan
310	index = imaginary(points)#replicate(1, nzjpt)+replicate(nxny, npoints)#lindgen(nz*jpt)
311	index = reform(index, npoints, nz, jpt, /over)
312	value2 = array[index]
313	terre = where(value2 GT valmask/10)
314	if terre[0] NE -1 then value2[terre] = !values.f_nan
315	poids = poids#replicate(1, nz*jpt)
316	poids = reform(poids, npoints, nz, jpt, /over)
317	res = poidsvalue1+(1-poids)value2
318	; moyenne suivant z ?
319	if strpos(type, 'z') EQ -1 then begin
320	nan = where(finite(res) EQ 0)
321	if vargrid EQ 'W' then e3 = e3w[firstzw:lastzw] ELSE e3 = e3t[firstzt:lastzt]
322	weight = replicate(1, npoints)#e3
323	weight = weight[*]#replicate(1, jpt)
324	weight = reform(weight, npoints, nz, jpt, /over)
325	if nan[0] NE -1 then weight[nan] = !values.f_nan
326	totalweight = total(weight, 2, /nan)
327	zero = where(totalweight EQ 0)
328	if zero[0] NE -1 then totalweight[zero] = !values.f_nan
329	res = total(res*weight, 2, /nan)/totalweight
330	direc = 'z'+string(byte(testvar(var = toto)))
331	endif
332	END
333	endcase
334	;---------------------
335
336	terre = where(finite(res) EQ 0)
337	if terre[0] NE -1 then res[terre] = valmask
338
339	if n_elements(showbuild) then BEGIN
340	winsave = !window
341	psave = !p
342	xsave = !x
343	ysave = !y
344	plt, findgen(nx, ny), /nodata, /nofill, /rempli, title = '', subtitle = '', coast_thick = 2, window = showbuild
345	!p.title = ''
346	!p.subtitle = ''
347
348	plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50
349	plots, [endpoints[0], endpoints[2]], [endpoints[1], endpoints[3]], color = 50, psym = 2, thick = 2
350
351	FOR i = 0, n_elements(points1)-1 DO $
352	plots, [float(points1[i]), float(points2[i])] $
353	, [imaginary(points1[i]), imaginary(points2[i])], color = 150
354
355	plots, float(points1), imaginary(points1), color = 150, psym = 1
356	plots, float(points2), imaginary(points2), color = 150, psym = 1
357	plots, float(inter), imaginary(inter), color = 250, psym = 1
358
359	IF terre[0] NE -1 THEN plots, float(inter[terre]), imaginary(inter[terre]), color = 0, psym = 1
360
361	; dummy = ''
362	; read, dummy, prompt = 'press return to continue'
363	IF !d.name EQ 'PS' THEN erase ELSE wset, winsave
364	!p = psave
365	!x = xsave
366	!y = ysave
367	ENDIF
368
369	restoreboxparam, 'boxparam4section.dat'
370
371	;---------------------
372	return
373	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format