- Timestamp:
- 05/02/06 15:32:01 (18 years ago)
- Location:
- trunk
- Files:
-
- 1 deleted
- 13 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/ToBeReviewed/TRIANGULATION/ciseauxtri.pro
r27 r29 38 38 ;------------------------------------------------------------ 39 39 FUNCTION ciseauxtri, triang, glam, gphi, TOUT = tout, _EXTRA = ex 40 @common 40 ;--------------------------------------------------------- 41 @cm_4mesh 42 IF NOT keyword_set(key_forgetold) THEN BEGIN 43 @updatenew 44 ENDIF 45 ;--------------------------------------------------------- 46 IF NOT keyword_set(key_periodic) AND NOT keyword_set(key_irregular) $ 47 AND NOT (!map.projection LE 7 AND !map.projection NE 0) $ 48 AND NOT (!map.projection EQ 14 OR !map.projection EQ 15 $ 49 OR !map.projection EQ 18) THEN return, triang 50 ; 41 51 tempsun = systime(1) ; pour key_performance 52 ; 42 53 taille = size(glam) 43 54 nx = taille[1] … … 45 56 46 57 tempdeux = systime(1) ; pour key_performance =2 47 z = convert_coord(glam (*),gphi(*),/data,/to_normal)48 x = z (0, *)49 y = z (1, *)58 z = convert_coord(glam[*],gphi[*],/data,/to_normal) 59 x = z[0, *] 60 y = z[1, *] 50 61 tempvar = SIZE(TEMPORARY(z)) ; delete z 51 62 IF testvar(var = key_performance) EQ 2 THEN $ … … 60 71 OR !map.projection EQ 14 OR !map.projection EQ 15 OR !map.projection EQ 18 then begin 61 72 tempdeux = systime(1) ; pour key_performance =2 62 63 ind = where(( finite((x*y)[triang[0, *]])$ ; points NaN ds la premiere colonne 64 *finite((x*y)[triang[1, *]]) $ ; points NaN ds la 2eme colonne 65 *finite((x*y)[triang[2, *]]) ) EQ 1) ; points NaN ds la 3eme colonne 66 67 if ind[0] NE -1 then triang = triang[*, ind] ELSE return, -1 73 ; 74 test = (x*y)[triang] 75 test = finite(temporary(test), /nan) 76 test = total(temporary(test), 1) 77 ind = where(temporary(test) EQ 0) 78 ; 79 if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return, -1 80 trichanged = 1b 68 81 IF testvar(var = key_performance) EQ 2 THEN $ 69 82 print, 'temps ciseauxtri: recherche points a NAN', systime(1)-tempdeux 70 83 endif 71 84 ; 72 73 seuil = 5. 85 seuil = 5 < (min([nx, ny])-2) 74 86 ; 75 87 ; maintenant on vire les triangles dont un des cotes a une taille … … 77 89 ; 78 90 79 if keyword_set(key_periodique) then begin 80 tempdeux = systime(1) ; pour key_performance =2 81 ind=where( (abs(x[triang[1,*]]-x[triang[0,*]]) LE (!p.position[2]-!p.position[0])/seuil) $ 82 AND (abs(x[triang[2,*]]-x[triang[1,*]]) LE (!p.position[2]-!p.position[0])/seuil) $ 83 AND (abs(x[triang[0,*]]-x[triang[2,*]]) LE (!p.position[2]-!p.position[0])/seuil) $ 84 ; 85 AND (abs(y[triang[0,*]]-y[triang[1,*]]) LE (!p.position[3]-!p.position[1])/seuil) $ 86 AND (abs(y[triang[1,*]]-y[triang[2,*]]) LE (!p.position[3]-!p.position[1])/seuil) $ 87 AND (abs(y[triang[2,*]]-y[triang[0,*]]) LE (!p.position[3]-!p.position[1])/seuil) ) 88 91 if keyword_set(key_periodic) then begin 92 tempdeux = systime(1) ; pour key_performance =2 93 ; 94 xtri = x[triang] 95 xtri = xtri-shift(xtri, 1, 0) 96 testx = abs(temporary(xtri)) GT ((!p.position[2]-!p.position[0])/seuil) 97 testx = total(temporary(testx), 1) 98 ; 99 ytri = y[triang] 100 ytri = ytri-shift(ytri, 1, 0) 101 testy = abs(temporary(ytri)) GT ((!p.position[3]-!p.position[1])/seuil) 102 testy = total(temporary(testy), 1) 103 ; 104 test = temporary(testx)+temporary(testy) 105 ind=where(temporary(test) EQ 0) 106 ; 89 107 IF testvar(var = key_performance) EQ 2 THEN $ 90 108 print, 'temps ciseauxtri: trouver les triangles trop grands', systime(1)-tempdeux 91 109 ; 92 110 tempdeux = systime(1) ; pour key_performance =2 93 if ind[0] NE -1 then triang = triang[*, ind] ELSE return, -1 111 if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return, -1 112 trichanged = 1b 94 113 IF testvar(var = key_performance) EQ 2 THEN $ 95 114 print, 'temps ciseauxtri: virer les triangles trop grands', systime(1)-tempdeux … … 99 118 ; valable qd /TOUT est active 100 119 ; 101 ;if keyword_set(key_irregular) then begin120 if keyword_set(key_irregular) then begin 102 121 103 ; tempdeux = systime(1) ; pour key_performance =2 104 ; seuil = 0 105 ; seuil = .1*max([!p.position[2]-!p.position[0],!p.position[3]-!p.position[1]]) 106 ; ind=where(x[triang[0,*]] GE !p.position[0]-seuil AND x[triang[0,*]] LE !p.position[2]+seuil $ 107 ; AND x[triang[1,*]] GE !p.position[0]-seuil AND x[triang[1,*]] LE !p.position[2]+seuil $ 108 ; AND x[triang[2,*]] GE !p.position[0]-seuil AND x[triang[2,*]] LE !p.position[2]+seuil $ 109 ; ; 110 ; AND y[triang[0,*]] GE !p.position[1]-seuil AND y[triang[0,*]] LE !p.position[3]+seuil $ 111 ; AND y[triang[1,*]] GE !p.position[1]-seuil AND y[triang[1,*]] LE !p.position[3]+seuil $ 112 ; AND y[triang[2,*]] GE !p.position[1]-seuil AND y[triang[2,*]] LE !p.position[3]+seuil ) 113 ; IF testvar(var = key_performance) EQ 2 THEN $ 114 ; print, 'temps ciseauxtri: trouver les triangles hors du cadre', systime(1)-tempdeux 115 ; ; 116 ; tempdeux = systime(1) ; pour key_performance =2 117 ; if ind[0] NE -1 then triang = triang[*, ind] ELSE return, -1 118 ; IF testvar(var = key_performance) EQ 2 THEN $ 119 ; print, 'temps ciseauxtri: virer les triangles hors du cadre', systime(1)-tempdeux 120 ; ENDIF 122 tempdeux = systime(1) ; pour key_performance =2 123 xtri = x[triang] 124 test1 = xtri GE !p.position[0] 125 test2 = xtri LE !p.position[2] 126 undefine, xtri 127 testx = temporary(test1)*temporary(test2) 128 testx = total(temporary(testx), 1) 129 ; 130 ytri = y[triang] 131 test1 = ytri GE !p.position[1] 132 test2 = ytri LE !p.position[3] 133 undefine, ytri 134 testy = temporary(test1)*temporary(test2) 135 testy = total(temporary(testy), 1) 136 ; 137 test = temporary(testx)*temporary(testy); 138 ; 139 ind=where(temporary(test) NE 0) 140 ; 141 if ind[0] NE -1 then triang = triang[*, temporary(ind)] ELSE return, -1 142 trichanged = 1b 143 IF testvar(var = key_performance) EQ 2 THEN $ 144 print, 'temps ciseauxtri: virer les triangles hors du cadre', systime(1)-tempdeux 145 ENDIF 121 146 ; 122 147 ; dernier tri 123 148 ; 124 if n_elements(ind) GT 1then BEGIN149 if keyword_set(trichanged) then BEGIN 125 150 ; 126 151 ; il faut verifier que les triangles que l''on garde ne … … 149 174 ; ont des indices suivant x qui sont 0 et nx-1. Ils appatienent au 150 175 ; rectangles de la colonne nx-1 et non de la colonne 0 151 if keyword_set(key_periodi que) AND nx EQ jpi then BEGIN176 if keyword_set(key_periodic) AND nx EQ jpi then BEGIN 152 177 indxmax = indxtriang0 > indxtriang2 153 178 indxtriang = indxmin + (nx-1)*(indxmin EQ 0 AND indxmax EQ (nx-1)) … … 199 224 ; 200 225 if keyword_set(key_performance) THEN print, 'temps ciseauxtri', systime(1)-tempsun 226 ; 201 227 return, triang 202 228 end -
trunk/ToBeReviewed/TRIANGULATION/completecointerre.pro
r27 r29 16 16 ; KEYWORD PARAMETERS: _EXTRA 17 17 ; 18 ; CONT_COLOR: the color of the continent. defaut value is 19 ; (!d.n_colors - 1) < 255 => white 20 ; 18 21 ; OUTPUTS: non 19 22 ; … … 32 35 ;------------------------------------------------------------ 33 36 ;------------------------------------------------------------ 34 PRO completecointerre, COINMONTE = coinmonte, COINDESCEND = coindescend, INDICEZOOM = indicezoom $ 37 PRO draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 38 @cm_4mesh 39 ; the triangle must not be out of the domain 40 IF min(lons, max = maxlon) GE lon1 AND maxlon LE lon2 $ 41 AND min(lats, max = maxlat) GE lat1 AND maxlat LE lat2 then BEGIN 42 ; the triangle must not be too big 43 z = convert_coord(lons, lats, /data, /to_normal) 44 alldist = [(z[0, 2]-z[0, 0])^2 + (z[1, 2]-z[1, 0])^2 $ 45 , (z[0, 0]-z[0, 1])^2 + (z[1, 0]-z[1, 1])^2 $ 46 , (z[0, 1]-z[0, 2])^2 + (z[1, 1]-z[1, 2])^2] 47 IF max(alldist) LT seuil^2 THEN polyfill, lons, lats $ 48 , color = cont_color, _extra = ex 49 return 50 ENDIF 51 end 52 ;------------------------------------------------------------ 53 ;------------------------------------------------------------ 54 PRO completecointerre, COINMONTE = coinmonte, COINDESCEND = coindescend $ 55 , CONT_COLOR = cont_color, INDICEZOOM = indicezoom $ 35 56 , _extra = ex 36 57 @common 37 58 ;------------------------------------------------------------ 38 tempsun = systime(1) ; pour key_performance 59 ; if NOT keyword_set(coinmonte) then return 60 ; if NOT keyword_set(coindescend) then return 61 ; if NOT keyword_set(indicezoom) then return 62 tempsun = systime(1) ; pour key_performance 39 63 ;------------------------------------------------------------ 40 64 ; definitions des vecteurs coinmont et coindesc 41 65 ;------------------------------------------------------------ 42 if keyword_set(coinmonte) then coinmont = coinmonte $ 43 ELSE coinmont = cointerremont 44 if keyword_set(coindescend) then coindesc = coindescend $ 45 ELSE coindesc = cointerredesc 66 if keyword_set(coinmonte) then coinmont = coinmonte $ 67 ELSE coinmont = twin_corners_up 68 if keyword_set(coindescend) then coindesc = coindescend $ 69 ELSE coindesc = twin_corners_dn 70 IF NOT keyword_set(cont_color) THEN cont_color = (!d.n_colors-1) < 255 46 71 ;------------------------------------------------------------ 47 72 ; definition descoordonnees des points numerotes 1,2,3,4,5,6 cf. les 48 73 ; schemas en dessous! 49 74 ;------------------------------------------------------------ 50 tempdeux = systime(1); pour key_performance =251 52 if keyword_set(indicezoom) then begin53 long1 = glamv[indicezoom] & lati1 = gphiv[indicezoom] 54 long2 = glamu[indicezoom] & lati2 = gphiu[indicezoom] 55 long3 = glamf[indicezoom] & lati3 = gphif[indicezoom] 56 taille = size(indicezoom)57 nx = taille[1]75 tempdeux = systime(1) ; pour key_performance =2 76 if coinmont[0] NE -1 OR coindesc[0] NE -1 then BEGIN 77 if keyword_set(indicezoom) then BEGIN 78 ; if we use key_stide, the t, u, v and f points are no more related to 79 ; the same cell because glamf and gphif has be recomputed to be in the 80 ; middle of two t points 81 IF total(key_stride) EQ 3 AND finite(glamv[0]*gphiv[0]) NE 0 THEN BEGIN 82 long1 = glamv[indicezoom] & lati1 = gphiv[indicezoom] 58 83 ENDIF ELSE BEGIN 59 long1 = glamv & lati1 = gphiv 60 long2 = glamu & lati2 = gphiu 61 long3 = glamf & lati3 = gphif 62 nx = jpi 63 ENDELSE 64 long4 = long2 & lati4 = lati2 65 long5 = long1 & lati5 = lati1 66 endif 67 IF testvar(var = key_performance) EQ 2 THEN $ 84 long1 = glamt[indicezoom] & lati1 = gphif[indicezoom] 85 ENDELSE 86 IF total(key_stride) EQ 3 AND finite(glamu[0]*gphiu[0]) NE 0 THEN BEGIN 87 long2 = glamu[indicezoom] & lati2 = gphiu[indicezoom] 88 ENDIF ELSE BEGIN 89 long2 = glamf[indicezoom] & lati2 = gphit[indicezoom] 90 ENDELSE 91 long3 = glamf[indicezoom] & lati3 = gphif[indicezoom] 92 ENDIF ELSE BEGIN 93 IF total(key_stride) EQ 3 AND finite(glamv[0]*gphiv[0]) NE 0 THEN BEGIN 94 long1 = glamv & lati1 = gphiv 95 ENDIF ELSE BEGIN 96 long1 = glamt & lati1 = gphif 97 ENDELSE 98 IF total(key_stride) EQ 3 AND finite(glamu[0]*gphiu[0]) NE 0 THEN BEGIN 99 long2 = glamu & lati2 = gphiu 100 ENDIF ELSE BEGIN 101 long2 = glamf & lati2 = gphit 102 ENDELSE 103 long3 = glamf & lati3 = gphif 104 ENDELSE 105 ; 106 nx = (size(long1, /dimensions))[0] 107 ny = (size(long1, /dimensions))[1] 108 seuil = 5 < (min([nx, ny])-2) 109 seuil = min([(!p.position[2]-!p.position[0])/seuil $ 110 , (!p.position[3]-!p.position[1])/seuil]) 111 ; 112 ENDIF 113 ; 114 IF testvar(var = key_performance) EQ 2 THEN $ 68 115 print, 'temps completecointerre: positions des points', systime(1)-tempdeux 69 116 ; … … 75 122 ; 4 76 123 ; t(i+nx)=1 u(i+nx) t(i+nx+1)=0 77 ; | 78 ; | 124 ; | \ 125 ; | \ 79 126 ; 1 3 | \ 5 80 127 ; v(i)---------f(i)------------v(i+1) … … 85 132 ; 86 133 ; 87 if coinmont[0] NE -1 then BEGIN 88 tempdeux = systime(1) ; pour key_performance =2 89 for id = 0, n_elements(coinmont)-1 do BEGIN 90 i = coinmont[id] 91 if long1[i] GE lon1 AND long5[i+1] LE lon2 $ 92 AND lati2[i] GE lat1 AND lati4[i+nx] LE lat2 then begin 93 polyfill, [long1[i], long2[i], long3[i], long4[i+nx], long5[i+1], long3[i]] $ 94 , [lati1[i], lati2[i], lati3[i], lati4[i+nx], lati5[i+1], lati3[i]] $ 95 , color = c_cont, _extra = ex 96 endif 97 endfor 98 IF testvar(var = key_performance) EQ 2 THEN $ 99 print, 'temps completecointerre: trace de cointerremonte', systime(1)-tempdeux 100 endif 134 if coinmont[0] NE -1 then BEGIN 135 tempdeux = systime(1) ; pour key_performance =2 136 for id = 0, n_elements(coinmont)-1 do BEGIN 137 i = coinmont[id] 138 ii = i MOD nx 139 ij = i/nx 140 ; bottom triangle 141 lons = [long1[i], long2[i], long3[i]] 142 lats = [lati1[i], lati2[i], lati3[i]] 143 draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 144 ; upper triangle 145 IF ii NE nx-1 AND ij NE ny-1 THEN BEGIN 146 lons = [long3[i], long1[i+1], long2[i+nx]] 147 lats = [lati3[i], lati1[i+1], lati2[i+nx]] 148 draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 149 ENDIF 150 ENDFOR 151 IF testvar(var = key_performance) EQ 2 THEN $ 152 print, 'temps completecointerre: trace de cointerremonte', systime(1)-tempdeux 153 ENDIF 101 154 ;------------------------------------------------------------ 102 155 ; cas coin terre en descendante.: … … 115 168 ; t(i)=0 2 u(i) t(i+1)=1 116 169 ; 117 if keyword_set(coindescend) then coindesc = coindescend $ 118 ELSE coindesc = cointerredesc 119 if coindesc[0] NE -1 then begin 120 tempdeux = systime(1) ; pour key_performance =2 121 for id = 0, n_elements(coindesc)-1 do BEGIN 122 i = coindesc[id] 123 if long1[i] GE lon1 AND long5[i+1] LE lon2 $ 124 AND lati2[i] GE lat1 AND lati4[i+nx] LE lat2 then begin 125 polyfill, [long1[i], long4[i+nx], long3[i], long2[i], long5[i+1], long3[i]] $ 126 , [lati1[i], lati4[i+nx], lati3[i], lati2[i], lati5[i+1], lati3[i]] $ 127 , color = c_cont, _extra = ex 128 endif 129 endfor 130 IF testvar(var = key_performance) EQ 2 THEN $ 131 print, 'temps completecointerre: trace de cointerredescend', systime(1)-tempdeux 132 endif 170 if coindesc[0] NE -1 then begin 171 tempdeux = systime(1) ; pour key_performance =2 172 for id = 0, n_elements(coindesc)-1 do BEGIN 173 i = coindesc[id] 174 ii = i MOD nx 175 ij = i/nx 176 IF ii NE nx-1 AND ij NE ny-1 THEN BEGIN 177 ; left triangle 178 lons = [long1[i], long3[i], long2[i+nx]] 179 lats = [lati1[i], lati3[i], lati2[i+nx]] 180 draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 181 ; right triangle 182 lons = [long3[i], long2[i], long1[i+1]] 183 lats = [lati3[i], lati2[i], lati1[i+1]] 184 draw_corner_triangle, lons, lats, seuil, CONT_COLOR = cont_color, _extra = ex 185 ENDIF 186 ENDFOR 187 IF testvar(var = key_performance) EQ 2 THEN $ 188 print, 'temps completecointerre: trace de cointerredescend', systime(1)-tempdeux 189 ENDIF 133 190 134 191 ;------------------------------------------------------------ 135 192 IF keyword_set(key_performance) THEN print, 'temps completecointerre', systime(1)-tempsun 136 193 ;------------------------------------------------------------ 137 194 return 138 195 end -
trunk/ToBeReviewed/TRIANGULATION/dessinetri.pro
r27 r29 22 22 ; triangulate) 23 23 ; 24 ; KEYWORD PARAMETERS: tous ceux de plots 24 ; KEYWORD PARAMETERS: 25 ; 26 ; All plots or polyfill keywords. 27 ; 28 ; WAIT=x. to call wait x second between each triangle draw. 29 ; 30 ; /ONEBYONE: to draw the triangles one by one 31 ; 32 ; /FILL: to fill the triangles (using polyfill) instead of plotting them 33 ; 34 ; CHANGECOLOR=n. to change the color of each traingle. n colors 35 ; will be used and repeted if necessary. 36 ; 25 37 ; 26 38 ; OUTPUTS: … … 41 53 ;------------------------------------------------------------ 42 54 43 PRO dessinetri, tri, x, y, _extra = ex55 PRO dessinetri, tri, x, y, WAIT = wait, ONEBYONE = onebyone, FILL = fill, CHANGECOLOR = changecolor, _extra = ex 44 56 @common 45 57 tempsun = systime(1) ; pour key_performance 46 47 if n_params() EQ 3 then begin 48 glam = x 49 gphi = y 58 a = '' 59 if n_params() EQ 3 then BEGIN 60 CASE size(x, /n_dimensions)+size(y, /n_dimensions) OF 61 2:BEGIN 62 nx = n_elements(x) 63 ny = n_elements(y) 64 glam = x#replicate(1., ny) 65 gphi = replicate(1., nx)#y 66 END 67 4:BEGIN 68 glam = x 69 gphi = y 70 END 71 ELSE:BEGIN 72 dummy = report('x and y inputs of dessinetri must have the same number of dimensions (1 or 2)') 73 return 74 END 75 ENDCASE 50 76 ENDIF ELSE BEGIN 51 77 grille,mask,glam,gphi, tri = tri … … 53 79 tri = ciseauxtri(tri, glam, gphi) 54 80 ENDELSE 55 81 ; 82 IF keyword_set(changecolor) THEN BEGIN 83 oldname = !d.name 84 if !d.name EQ 'PS' OR !d.name EQ 'Z' then BEGIN 85 thisos = strupcase(strmid(!version.os_family, 0, 3)) 86 CASE thisOS of 87 'MAC': set_plot, thisOS 88 'WIN': set_plot, thisOS 89 ELSE: set_plot, 'X' 90 ENDCASE 91 ncolors=(!d.n_colors-1) < 255 92 set_plot, oldname 93 ENDIF ELSE ncolors=(!d.n_colors-1) < 255 94 color = 1+indgen(changecolor)*(ncolors/(changecolor-1)) 95 ENDIF ELSE color = 0 96 color = color#replicate(1, n_elements(tri)/3/n_elements(color)+1) 97 ; 56 98 tempdeux = systime(1) ; pour key_performance =2 57 99 for i = 0L, n_elements(tri)/3-1 do begin 58 100 t = [tri[*, i], tri[0, i]] 59 ; plots, glam[t], gphi[t], color = i MOD 255, _extra = ex 60 ; wait, .1 61 plots, glam[t], gphi[t], color = 0, _extra = ex 101 IF keyword_set(fill) THEN $ 102 polyfill, glam[t], gphi[t], color = color[i], _extra = ex $ 103 ELSE plots, glam[t], gphi[t], color = color[i], _extra = ex 104 IF keyword_set(wait) THEN wait, wait 105 IF keyword_set(onebyone) THEN read, a, prompt = 'press a key' 62 106 ENDFOR 63 107 IF testvar(var = key_performance) EQ 2 THEN $ -
trunk/ToBeReviewed/TRIANGULATION/drawcoast_c.pro
r27 r29 1 PRO drawcoast_c, mask, xf, yf, nx, ny, CONT_THICK = cont_thick, YSEUIL = yseuil, XSEUIL = xseuil, _extra = ex 2 @common 1 PRO drawcoast_c, mask, xf, yf, nx, ny, COAST_COLOR = coast_color, COAST_THICK = coast_thick, YSEUIL = yseuil, XSEUIL = xseuil, _extra = ex 2 ;--------------------------------------------------------- 3 @cm_4mesh 4 IF NOT keyword_set(key_forgetold) THEN BEGIN 5 @updatenew 6 @updatekwd 7 ENDIF 8 ;--------------------------------------------------------- 3 9 tempsun = systime(1) ; pour key_performance 4 10 ;--------------------------------------------------------- … … 6 12 ; on trace les segments verticaux: 7 13 ; 8 if NOT keyword_set( xseuil) then xseuil = 5.9 distanceseuil = (!p.position[ 2]-!p.position[0])/xseuil14 if NOT keyword_set(yseuil) then yseuil = 5. < (min([nx, ny])-2) 15 distanceseuil = (!p.position[3]-!p.position[1])/yseuil 10 16 ; liste: liste des points i pourlesquels on va tracer un segment entre 11 17 ; le point i,j-1 et i,j 12 18 tempdeux = systime(1) ; pour key_performance =2 13 liste = where( mask+shift(mask, -1, 0) EQ 1 $14 AND ( xf-shift(xf, 0, 1))^2+(yf-shift(yf, 0, 1))^2 LE distanceseuil)19 liste = where((mask+shift(mask, -1, 0)) EQ 1 $ 20 AND ((xf-shift(xf, 0, 1))^2+(yf-shift(yf, 0, 1))^2) LE distanceseuil^2) 15 21 IF liste[0] NE -1 THEN BEGIN 16 22 ; on recupere lx et ly qui sont les indices ds un tableau 2d des … … 26 32 print, 'temps tracecote: determiner liste des points concernes par un trait vertical', systime(1)-tempdeux 27 33 tempdeux = systime(1) ; pour key_performance =2 28 for pt = 0 , n_elements(lx)-1 do BEGIN34 for pt = 0L, n_elements(lx)-1 do BEGIN 29 35 i = lx[pt] & j = ly[pt] 30 36 plots, [xf[i, j-1], xf[i, j]], [yf[i, j-1], yf[i, j]] $ 31 , color=c_cote,thick=cont_thick, /normal, _extra = ex37 , color = coast_color, thick = coast_thick, /normal, _extra = ex 32 38 endfor 33 39 IF testvar(var = key_performance) EQ 2 THEN $ … … 43 49 ; periodique mais pour le plots 44 50 tempdeux = systime(1) ; pour key_performance =2 45 if keyword_set(key_periodi que) AND nx EQ jpi then begin51 if keyword_set(key_periodic) AND nx EQ jpi then begin 46 52 mask = [mask, mask[0, *]] 47 53 xf = [xf, xf[0, *]] … … 49 55 nx = nx+1 50 56 ENDIF 51 if NOT keyword_set( yseuil) then yseuil = 5.52 distanceseuil = (!p.position[ 3]-!p.position[1])/yseuil53 liste = where( mask+shift(mask, 0, -1) EQ 1 $54 AND ( xf-shift(xf, 1, 0))^2+(yf-shift(yf, 1, 0))^2 LE distanceseuil)57 if NOT keyword_set(xseuil) then xseuil = 5. < (min([nx, ny])-2) 58 distanceseuil = (!p.position[2]-!p.position[0])/xseuil 59 liste = where((mask+shift(mask, 0, -1)) EQ 1 $ 60 AND ((xf-shift(xf, 1, 0))^2+(yf-shift(yf, 1, 0))^2) LE distanceseuil^2) 55 61 IF liste[0] NE -1 THEN BEGIN 56 62 ly = liste/nx & lx = temporary(liste)-nx*ly … … 63 69 print, 'temps tracecote: determiner liste des points concernes par un trait horizontal', systime(1)-tempdeux 64 70 tempdeux = systime(1) ; pour key_performance =2 65 for pt = 0 , n_elements(lx)-1 do BEGIN71 for pt = 0L, n_elements(lx)-1 do BEGIN 66 72 i = lx[pt] & j = ly[pt] 67 73 plots, [xf[i-1, j], xf[i, j]], [yf[i-1, j], yf[i, j]] $ 68 , color=c_cote,thick=cont_thick, /normal, _extra = ex74 , color = coast_color, thick = coast_thick, /normal, _extra = ex 69 75 endfor 70 76 IF testvar(var = key_performance) EQ 2 THEN $ -
trunk/ToBeReviewed/TRIANGULATION/drawcoast_e.pro
r27 r29 1 PRO drawcoast_e, mask, xf, yf, nx, ny, CONT_THICK = cont_thick, YSEUIL = yseuil, XSEUIL = xseuil, onemore = onemore, _extra = ex 2 @common 1 PRO drawcoast_e, mask, xf, yf, nx, ny, COAST_COLOR = coast_color, COAST_THICK = coast_thick, YSEUIL = yseuil, XSEUIL = xseuil, onemore = onemore, _extra = ex 2 ;--------------------------------------------------------- 3 @cm_4mesh 4 IF NOT keyword_set(key_forgetold) THEN BEGIN 5 @updatenew 6 @updatekwd 7 ENDIF 8 ;--------------------------------------------------------- 3 9 tempsun = systime(1) ; pour key_performance 4 10 ;--------------------------------------------------------- 5 if keyword_set(key_periodi que) AND nx EQ jpi then begin11 if keyword_set(key_periodic) AND nx EQ jpi then begin 6 12 mask = [mask, mask[0, *]] 7 13 xf = [xf, xf[0, *]] … … 13 19 ; 14 20 if NOT keyword_set(onemore) then onemore = 0 15 if NOT keyword_set(xseuil) then xseuil = 5. 21 if NOT keyword_set(xseuil) then xseuil = 5. < (min([nx, ny])-2) 16 22 distanceseuil = (!p.position[2]-!p.position[0])/xseuil 17 23 ; liste: liste des points i pourlesquels on va tracer un segment … … 20 26 indexbis = index-nx+((index/nx+onemore) MOD 2) 21 27 liste = where(mask[index+1]+mask[indexbis] EQ 1 $ 22 AND (xf[index]-xf[indexbis])^2+(yf[index]-yf[indexbis])^2 LE distanceseuil )28 AND (xf[index]-xf[indexbis])^2+(yf[index]-yf[indexbis])^2 LE distanceseuil^2) 23 29 IF liste[0] NE -1 THEN BEGIN 24 30 index = index[liste] … … 26 32 for pt = 0, n_elements(index)-1 do begin 27 33 plots, [xf[index[pt]], xf[indexbis[pt]]], [yf[index[pt]], yf[indexbis[pt]]] $ 28 , color=c_cote,thick=cont_thick, /normal, _extra = ex34 , color = coast_color, thick = coast_thick, /normal, _extra = ex 29 35 endfor 30 36 ENDIF … … 32 38 ; we plot the borders of the diamond in this sense : / 33 39 ; 34 if NOT keyword_set(xseuil) then xseuil = 5. 40 if NOT keyword_set(xseuil) then xseuil = 5. < (min([nx, ny])-2) 35 41 distanceseuil = (!p.position[2]-!p.position[0])/xseuil 36 42 ; liste: liste des points i pourlesquels on va tracer un segment … … 39 45 indexbis = index+nx+((index/nx+onemore) MOD 2) 40 46 liste = where(mask[index+1]+mask[indexbis] EQ 1 $ 41 AND (xf[index]-xf[indexbis])^2+(yf[index]-yf[indexbis])^2 LE distanceseuil )47 AND (xf[index]-xf[indexbis])^2+(yf[index]-yf[indexbis])^2 LE distanceseuil^2) 42 48 IF liste[0] NE -1 THEN BEGIN 43 49 index = index[liste] … … 45 51 for pt = 0, n_elements(index)-1 do begin 46 52 plots, [xf[index[pt]], xf[indexbis[pt]]], [yf[index[pt]], yf[indexbis[pt]]] $ 47 , color=c_cote,thick=cont_thick, /normal, _extra = ex53 , color = coast_color, thick = coast_thick, /normal, _extra = ex 48 54 endfor 49 55 ENDIF -
trunk/ToBeReviewed/TRIANGULATION/section.pro
r27 r29 32 32 ;------------------------------------------------------------ 33 33 34 PRO section, field, res, glamaxe, gphiaxe, ENDPOINTS = endpoints, BOITE = boite, TYPE = type, DIREC = direc 35 @common 36 ; 37 ;--------------------- 38 array = litchamp(field) 39 ;--------------------- 40 ; definition de boite en fonction de endpoints 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 41 52 ; puis redefinition du domaine 42 boite2d = [min([endpoints[0], endpoints[2]]), max([endpoints[0], endpoints[2]]), min([endpoints[1], endpoints[3]]), max([endpoints[1], endpoints[3]])] 43 ; 44 minprof = 0 45 profdefault = 200 46 if n_elements(type) EQ 0 then type = 'nothing' 47 Case N_Elements(Boite) OF 48 1:boite=[boite2d, minprof, boite[0]] 49 2:boite=[boite2d, boite[0],boite[1]] 50 Else:$ 51 if strpos(type, 'z') NE -1 THEN boite=[boite2d, minprof, profdefault] $ 52 ELSE boite=[boite2d, prof1, prof2] 53 ENDCASE 54 if strpos(type, 'z') NE -1 THEN BEGIN 55 !y.range = [boite[5], boite[4]] ;on garde les yranges (axe z) avant de changer la boite. 56 profmax = boite[5] 57 vraiboite = boite 58 if vargrid EQ 'W' then gdep = gdepw ELSE gdep = gdept 59 ; check with vertical grid limits (nearest level) 60 gwork = gdep 61 ; check the increse or decrese of the z axis 62 IF gwork[1] LE gwork[0] THEN gwork = reverse(gdep, 1) 63 niveauprof = where(gwork ge boite[5]) & niveauprof = niveauprof[0] 64 if niveauprof NE -1 then boite[5] = gwork[niveauprof]+1 65 ENDIF 66 domdef, boite, GRILLE=[vargrid], /findalways, _extra = ex 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 67 82 ; on redefinit lon1, ... au cas ou findalways ait ete utilise ds domdef 68 lon1 = min([endpoints[0], endpoints[2]]) 69 lon2 = max([endpoints[0], endpoints[2]]) 70 lat1 = min([endpoints[1], endpoints[3]]) 71 lat2 = max([endpoints[1], endpoints[3]]) 72 ; on recupere la grille sur la boite 73 grille, mask, glam, gphi, gdep, nx, ny,nz,premierx,premiery,premierz,dernierx, derniery, dernierz 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 74 124 ;--------------------- 75 125 ; on definit la triangulation qui va nous permetre de determiner la … … 77 127 ; aussi bien que sur la mer. suivant le sens de la section -plutot 78 128 ; longitude ou plutot latitude- on definit la facon de trianguler 79 if strpos(type, 'x') NE -1 then BEGIN 80 downward = (lindgen(nx, ny))[0:nx-2, 0:ny-2] 81 tri = definetri(nx, ny, (downward)[*]) 82 ENDIF ELSE tri = definetri(nx, ny) 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 83 144 ;--------------------- 84 145 ; equation de la droite suivant laquelle on fait la section 85 146 ;--------------------- 86 87 88 147 abc = linearequation(endpoints[0:1], endpoints[2:3]) 148 glamtri = glam[tri] 149 gphitri = gphi[tri] 89 150 ; quels sont les points de la triangulation qui sont au dessus et au 90 151 ; dessous de la droite ? 91 152 if abc[1] NE 0 THEN $ 92 153 test = temporary(gphitri) GE (-abc[0]/abc[1]*temporary(glamtri)-abc[2]/abc[1]) $ 93 94 95 154 ELSE test = temporary(glamtri) GE (-abc[1]/abc[0]*temporary(gphitri)-abc[2]/abc[0]) 155 156 zero123 = total(test, 1) 96 157 ; to keep: triangles de la triangulation qui sont a cheval sur la droite. 97 98 99 100 101 102 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] 103 164 ; quel est le sommet du triangle qui est seul d''un cote de la droite? 104 105 106 107 108 109 110 111 112 113 114 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)] 115 176 ; points1 le point du triangles qui est seul d''un cote de la droite. 116 177 ; point2 l''autre point du triangle de l''autre cote de la droite 117 118 119 120 121 122 123 124 125 126 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)] 127 188 ; points : complexe contenant les couples de points de part et 128 189 ; d''autre de la droite. Ils faut supprimer les doublons. 129 130 131 132 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)] 133 194 ; points1 les coordnnees du point du triangles qui est seul d''un cote de la droite. 134 195 ; point2 les coordnnees de l''autre point du triangle de l''autre cote de la droite 135 136 196 points1 = complex(glam[ double(points)], gphi[ double(points)]) 197 points2 = complex(glam[imaginary(points)], gphi[imaginary(points)]) 137 198 ; droites les equations des droites dont on cherche l''intersection 138 199 ; avec la section. 139 140 200 droites = linearequation(points1, points2) 201 inter = lineintersection(droites, abc#replicate(1, n_elements(points1) ) ) 141 202 142 203 ; les ccordonnes geographiques des points que l''on cherche sur la section. 143 144 204 glamaxe = float(inter) 205 gphiaxe = imaginary(inter) 145 206 ; on les range ds l''ordre croissant entre les bornes de la section 146 if strpos(type, 'x') NE -1 then BEGIN 147 sort = sort(glamaxe) 148 glamaxe = glamaxe[sort] 149 inbox = where(glamaxe GE lon1 AND glamaxe LE lon2) 150 glamaxe = glamaxe[inbox] 151 sort = sort[inbox] 152 gphiaxe = gphiaxe[sort] 153 ENDIF ELSE BEGIN 154 sort = sort(gphiaxe) 155 gphiaxe = gphiaxe[sort] 156 inbox = where(gphiaxe GE lat1 AND gphiaxe LE lat2) 157 gphiaxe = gphiaxe[inbox] 158 sort = sort[inbox] 159 glamaxe = glamaxe[sort] 160 ENDELSE 161 points = points[sort] 162 points1 = points1[sort] 163 points2 = points2[sort] 164 inter = inter[sort] 165 poids = abs(points2-inter)/abs(points2-points1) 166 167 ;--------------------- 168 array = fitintobox(array) 169 if array[0] EQ -1 THEN BEGIN 170 res = -1 171 return 172 ENDIF 173 if n_elements(valmask) EQ 0 THEN valmask = 1e20 174 taille=size(array) 175 if jpt GT 1 AND taille[0] GE 3 AND strpos(type, 't') EQ -1 then BEGIN 176 direc = 't' 177 array = grossemoyenne(array, 't') 178 taille=size(array) 179 jpt = 1 180 ENDIF 181 case 1 of 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 182 244 ;---------------------------------------------------------------------------- 183 245 ;xy 184 246 ;---------------------------------------------------------------------------- 185 186 187 188 189 190 191 192 193 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 = poids*value1+(1-poids)*value2 255 END 194 256 ;---------------------------------------------------------------------------- 195 257 ;xyz 196 258 ;---------------------------------------------------------------------------- 197 198 199 200 201 202 203 204 205 206 207 208 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 = poids*value1+(1-poids)*value2 209 271 ; moyenne suivant z ? 210 211 212 if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt]213 214 215 216 217 218 219 direc = 'z'+string(byte(testvar(var=toto)))220 221 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 222 284 ;---------------------------------------------------------------------------- 223 285 ;xyt 224 286 ;---------------------------------------------------------------------------- 225 226 227 228 229 230 231 232 233 234 235 236 237 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 = poids*value1+(1-poids)*value2 299 END 238 300 ;---------------------------------------------------------------------------- 239 301 ;xyzt 240 302 ;---------------------------------------------------------------------------- 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 303 taille[0] EQ 4:BEGIN 304 npoints = n_elements(points) 305 index = double(points)#replicate(1, nz*jpt)+replicate(nx*ny, 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, nz*jpt)+replicate(nx*ny, 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 = poids*value1+(1-poids)*value2 256 318 ; moyenne suivant z ? 257 if strpos(type, 'z') EQ -1 then begin 258 nan = where(finite(res) EQ 0) 259 if vargrid EQ 'W' then e3 = e3w[premierzw:dernierzw] ELSE e3 = e3t[premierzt:dernierzt] 260 weight = replicate(1, npoints)#e3 261 weight = weight[*]#replicate(1, jpt) 262 weight = reform(weight, npoints, nz, jpt, /over) 263 if nan[0] NE -1 then weight[nan] = !values.f_nan 264 totalweight = total(weight, 2, /nan) 265 zero = where(totalweight EQ 0) 266 if zero[0] NE -1 then totalweight[zero] = !values.f_nan 267 res = total(res*weight, 2, /nan)/totalweight 268 direc = 'z'+string(byte(testvar(var=toto))) 269 endif 270 END 271 endcase 272 ;--------------------- 273 274 terre = where(finite(res) EQ 0) 275 if terre[0] NE -1 then res[terre] = valmask 276 277 278 ; plt, tmask[*,*,0],/nocouleur, /rempli, title = '', subtitle = '' 279 ; !p.title='' 280 ; !p.subtitle='' 281 ; dessinetri 282 ; plot,[endpoints[0],endpoints[2]],[endpoints[1],endpoints[3]],/noerase,color=50 283 284 ; plot,float(points1),imaginary(points1),/noerase,color=150,psym=1 285 ; plot,float(points2),imaginary(points2),/noerase,color=150,psym=1 286 ; plot,float(inter),imaginary(inter),/noerase,color=250,psym=1 287 288 ; plot,[float(points1[0])],[imaginary(points1[0])],/noerase,color=150,psym=1 289 ; plot,[float(points2[0])],[imaginary(points2[0])],/noerase,color=150,psym=1 290 ; plot,[float(inter[0])],[imaginary(inter[0])],/noerase,color=250,psym=1 291 292 293 ;--------------------- 294 return 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 295 373 end -
trunk/ToBeReviewed/TRIANGULATION/tracecote.pro
r27 r29 15 15 ; KEYWORD PARAMETERS: 16 16 ; 17 ; CONT_THICK: l''epaisseur du trait pour tracer les 17 ; COAST_COLOR: the color of the coastline. 18 ; defaut value is 0 => black 19 ; 20 ; COAST_THICK: l''epaisseur du trait pour tracer les 18 21 ; continents. par defaut c''est 1. 22 ; 23 ; /SURFACE_COASTLINE: to draw the furface coast line instead of 24 ; the coast line at level firstz[tw]. Usefull only for deep 25 ; plots! 19 26 ; 20 27 ; XSEUIL: pour eliminer les segments de cote qui sont trop … … 45 52 ;------------------------------------------------------------ 46 53 ;------------------------------------------------------------ 47 ; PRO tracecote, CONT_THICK = cont_thick, YSEUIL = yseuil, XSEUIL = xseuil, _extra = ex 48 PRO tracecote, _extra = ex 49 @common 54 PRO tracecote, SURFACE_COASTLINE = surface_coastline, _EXTRA = ex 55 ;-------------------------------------------------------------- 56 ; include commons 57 @cm_4data 58 @cm_4mesh 59 IF NOT keyword_set(key_forgetold) THEN BEGIN 60 @updatenew 61 ENDIF 62 ;-------------------------------------------------------------- 50 63 tempsun = systime(1) ; pour key_performance 51 64 if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' … … 56 69 ; 57 70 tempdeux = systime(1) ; pour key_performance =2 58 premierx = 0 > (min([premierxt, premierxu])-1)59 dernierx = (max([dernierxt, dernierxu])+1) < (jpi-1)60 premiery = 0 > (min([premieryt, premieryv])-1)61 derniery = (max([dernieryt, dernieryv])+1) < (jpj-1)62 nx = dernierx-premierx+163 ny = derniery-premiery+171 firstx = 0 > (min([firstxt, firstxf])-1) 72 lastx = (max([lastxt, lastxf])+1) < (jpi-1) 73 firsty = 0 > (min([firstyt, firstyf])-1) 74 lasty = (max([lastyt, lastyf])+1) < (jpj-1) 75 nx = lastx-firstx+1 76 ny = lasty-firsty+1 64 77 ; quel niveau vertical choisir ? 65 IF strupcase(vargrid) eq 'W' THEN BEGIN 66 nz = nzw & premierz = premierzw & dernierz = dernierzw 67 ENDIF ELSE BEGIN 68 nz = nzt & premierz = premierzt & dernierz = dernierzt 69 ENDELSE 70 if nz eq jpk then nivz = niveau-1 ELSE nivz = nz-1 78 IF keyword_set(surface_coastline) THEN firstz = 0 ELSE $ 79 IF strupcase(vargrid) eq 'W' THEN firstz = firstzw ELSE firstz = firstzt 71 80 ; attribution du masque et des coordonnes delimitant les limites de la 72 81 ; terre (coordonnees f) 73 mask = tmask[ premierx:dernierx, premiery:derniery, niveau-1]74 xf = glamf[ premierx:dernierx, premiery:derniery]75 yf = gphif[ premierx:dernierx, premiery:derniery] ;82 mask = tmask[firstx:lastx, firsty:lasty, firstz] 83 xf = glamf[firstx:lastx, firsty:lasty] 84 yf = gphif[firstx:lastx, firsty:lasty] ; 76 85 IF testvar(var = key_performance) EQ 2 THEN $ 77 86 print, 'temps tracecote: determiner mask xf yf', systime(1)-tempdeux … … 107 116 tempvar = SIZE(TEMPORARY(ind)) ; on efface ind 108 117 ; 109 ; pour eviter d'avoir des traits de cotes qui travessent tout l''ecran110 ; (qui peuvent etre tres proches sur la sphere mais tres eloignes sur111 ; le dessin suivant la projection) on va selectionner uniquements les112 ; segments de trait de cote qui ne sont pas trop grand, cad qui ne113 ; depassent pas la taille de la fenetre/seuil et pour lesquel le114 ; maskque change de valeur entre les extremites du segment.115 ;116 ; on trace les segments verticaux:117 ;118 ; if NOT keyword_set(xseuil) then xseuil = 5.119 ; distanceseuil = (!p.position[2]-!p.position[0])/xseuil120 ; ; liste: liste des points i pourlesquels on va tracer un segment entre121 ; ; le point i,j-1 et i,j122 ; tempdeux = systime(1) ; pour key_performance =2123 ; liste = where(mask+shift(mask, -1, 0) EQ 1 $124 ; AND (xf-shift(xf, 0, 1))^2+(yf-shift(yf, 0, 1))^2 LE distanceseuil)125 ; IF liste[0] NE -1 THEN BEGIN126 ; ; on recupere lx et ly qui sont les indices ds un tableau 2d des127 ; ; points donnes par liste128 ; ly = liste/nx & lx = temporary(liste)-nx*ly129 ; indice = where(ly NE 0) ; on ne prend pas les points concernant130 ; if indice[0] NE -1 then begin131 ; ; la premiere ligne car ds ce cas le pt j-1 n''est pas definit132 ; lx = lx[indice] & ly = ly[temporary(indice)]133 ; ; boucle sur les points concernes et trace du segment134 ; ; rq: on utilise plost au lieu de plot car plots est bcp plus rapide.135 ; IF testvar(var = key_performance) EQ 2 THEN $136 ; print, 'temps tracecote: determiner liste des points concernes par un trait vertical', systime(1)-tempdeux137 ; tempdeux = systime(1) ; pour key_performance =2138 ; for pt = 0, n_elements(lx)-1 do BEGIN139 ; i = lx[pt] & j = ly[pt]140 ; plots, [xf[i, j-1], xf[i, j]], [yf[i, j-1], yf[i, j]] $141 ; , color=c_cote,thick=cont_thick, /normal, _extra = ex142 ; endfor143 ; IF testvar(var = key_performance) EQ 2 THEN $144 ; print, 'temps tracecote: trace des traits verticaux', systime(1)-tempdeux145 ; endif146 ; ENDIF147 ; ;148 ; ; pour le trace des segments horizontaux, c''est la meme chose sauf149 ; ; qu'il faut faire attention si on est periodique:150 ; ;151 ; ; si on est periodique on duplique la premiere colonne et on la met a152 ; ; la fin. (ceci est fait non pas pour le shift qui est par defaut153 ; ; periodique mais pour le plots154 ; tempdeux = systime(1) ; pour key_performance =2155 ; if keyword_set(key_periodique) AND nx EQ jpi then begin156 ; mask = [mask, mask[0, *]]157 ; xf = [xf, xf[0, *]]158 ; yf = [yf, yf[0, *]]159 ; nx = nx+1160 ; ENDIF161 ; if NOT keyword_set(yseuil) then yseuil = 5.162 ; distanceseuil = (!p.position[3]-!p.position[1])/yseuil163 ; liste = where(mask+shift(mask, 0, -1) EQ 1 $164 ; AND (xf-shift(xf, 1, 0))^2+(yf-shift(yf, 1, 0))^2 LE distanceseuil)165 ; IF liste[0] NE -1 THEN BEGIN166 ; ly = liste/nx & lx = temporary(liste)-nx*ly167 ; indice = where(ly NE ny-1 AND lx NE 0)168 ; if indice[0] NE -1 then begin169 ; ; on ne prend pas les points de la170 ; ; premiere colonne et de la derniere ligne (car on l''a rajoute artificiellement!))171 ; lx = lx[indice] & ly = ly[temporary(indice)]172 ; IF testvar(var = key_performance) EQ 2 THEN $173 ; print, 'temps tracecote: determiner liste des points concernes par un trait horizontal', systime(1)-tempdeux174 ; tempdeux = systime(1) ; pour key_performance =2175 ; for pt = 0, n_elements(lx)-1 do BEGIN176 ; i = lx[pt] & j = ly[pt]177 ; plots, [xf[i-1, j], xf[i, j]], [yf[i-1, j], yf[i, j]] $178 ; , color=c_cote,thick=cont_thick, /normal, _extra = ex179 ; endfor180 ; IF testvar(var = key_performance) EQ 2 THEN $181 ; print, 'temps tracecote: trace des traits horizontaux', systime(1)-tempdeux182 ; endif183 ; endif184 118 if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 185 119 case key_gridtype of -
trunk/ToBeReviewed/TRIANGULATION/tracemask.pro
r27 r29 17 17 ; KEYWORD PARAMETERS: 18 18 ; 19 ; CONT_THICK: l''epaisseur du trait pour tracer les 19 ; COAST_COLOR: the color of the coastline. 20 ; defaut value is 0 => black 21 ; 22 ; COAST_THICK: l''epaisseur du trait pour tracer les 20 23 ; continents. par defaut c''est 1. 21 24 ; … … 36 39 ;------------------------------------------------------------ 37 40 ;------------------------------------------------------------ 38 PRO tracemask, maskentree, xin, yin, CONT_THICK = cont_thick, OVERPLOT = overplot, _extra = ex 39 xentree = xin 40 yentree = yin 41 PRO tracemask, maskentree, xin, yin, COAST_COLOR = coast_color, COAST_THICK = coast_thick, OVERPLOT = overplot, _extra = ex 42 ; 41 43 if keyword_set(overplot) then return 42 @common 44 ;--------------------------------------------------------- 45 @cm_general 46 IF NOT keyword_set(key_forgetold) THEN BEGIN 47 @updatekwd 48 ENDIF 49 ;--------------------------------------------------------- 43 50 tempsun = systime(1) ; pour key_performance 44 51 ; on s''afranchit des problemes de bord: … … 47 54 nx = tailleentree[1]+1 48 55 ny = tailleentree[2]+1 49 ; on agrandi le mask de une colonne a gauche et de une colonne a droite. 56 ; we check the input axis 57 IF n_elements(xin) EQ 0 THEN xentree = findgen(nx-1) ELSE xentree = xin 58 IF (size(xentree))[0] EQ 1 THEN xentree = xentree#replicate(1,ny-1) 59 IF n_elements(yin) EQ 0 THEN yentree = findgen(ny-1) ELSE yentree = yin 60 IF (size(yentree))[0] EQ 1 THEN yentree = replicate(1,nx-1)#yentree 61 ; on agrandi le mask de une colonne a gauche et de une colonne en bas 50 62 mask = intarr(tailleentree[1]+1, tailleentree[2]+1) 51 63 mask[1:tailleentree[1], 1:tailleentree[2]] = maskentree … … 105 117 ; boucle sur les points concernes et trace du segment 106 118 ; rq: on utilise plots au lieu de plot car plots est bcp plus rapide. 107 for pt = 0 , n_elements(lx)-1 do BEGIN119 for pt = 0L, n_elements(lx)-1 do BEGIN 108 120 i = lx[pt] & j = ly[pt] 109 121 plots, [xf[i, j-1], xf[i, j]], [yf[i, j-1], yf[i, j]] $ 110 , color=c_cote,thick=cont_thick, _extra = ex122 , color = coast_color, thick = coast_thick, _extra = ex 111 123 if pt LT 5 then begin 112 124 endif … … 129 141 print, 'temps tracemask: liste traits horizontaux', systime(1)-tempdeux 130 142 tempdeux = systime(1) ; pour key_performance =2 131 for pt = 0 , n_elements(lx)-1 do BEGIN143 for pt = 0L, n_elements(lx)-1 do BEGIN 132 144 i = lx[pt] & j = ly[pt] 133 145 plots, [xf[i-1, j], xf[i, j]], [yf[i-1, j], yf[i, j]] $ 134 , color=c_cote,thick=cont_thick, _extra = ex146 , color = coast_color, thick = coast_thick, _extra = ex 135 147 endfor 136 148 IF testvar(var = key_performance) EQ 2 THEN $ -
trunk/ToBeReviewed/TRIANGULATION/triangule.pro
r27 r29 1 FUNCTION triangule, maskentree, REGULIER = regulier, _extra = ex1 FUNCTION triangule, maskentree, BASIC = basic, COINMONTE = coinmonte, COINDESCEND = coindescend, _extra = ex 2 2 @common 3 if keyword_set(regulier) then return, triangule_c(maskentree, REGULIER = regulier, _extra = ex) 3 ; 4 IF jpi EQ 1 OR jpj EQ 1 THEN return, -1 5 ; 6 IF arg_present(coinmonte) THEN coinmonte = 1 7 IF arg_present(coindescend) THEN coindescend = 1 8 ; 9 if keyword_set(basic) then $ 10 return, triangule_c(maskentree, /BASIC, COINMONTE = coinmonte $ 11 , COINDESCEND = coindescend, _extra = ex) 12 ; 4 13 if n_elements(key_gridtype) EQ 0 then key_gridtype = 'c' 5 14 if n_elements(maskentree) EQ 0 then maskentree = tmask[*, *, 0] 6 15 case key_gridtype of 7 16 'e':res = triangule_e(maskentree, _extra = ex) 8 'c':res = triangule_c(maskentree, REGULIER = regulier, _extra = ex)17 'c':res = triangule_c(maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend, _extra = ex) 9 18 endcase 10 19 return, res -
trunk/ToBeReviewed/TRIANGULATION/triangule_c.pro
r27 r29 42 42 ; 43 43 ; KEYWORD PARAMETERS: 44 ; /REGULIER: specifie que le masque est sur une grille reguliere 44 ; 45 ; /BASIC: specifie que le masque est sur une grille basice 45 46 ; (utiliser pour la triangulation ds les coupes verticales et 46 47 ; des hovmoellers) 48 ; 49 ; /KEEP_CONT: to keep the triangulation even on the continents 47 50 ; 48 51 ; COINMONTE=tableau, pour obtenir le tableau de "coins de terre 49 52 ; montant" a traiter avec completecointerre.pro ds la variable 50 53 ; tableau plutot que de la faire passer par la variable globale 51 ; cointerremont.54 ; twin_corners_up. 52 55 ; 53 56 ; COINDESCEND=tableau cf COINMONTE … … 78 81 ;------------------------------------------------------------ 79 82 ;------------------------------------------------------------ 80 FUNCTION triangule_c, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend $ 81 , REGULIER = regulier, PERIODIQUE = periodique 83 FUNCTION triangule_c, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend, BASIC = basic, KEEP_CONT = keep_cont 82 84 tempsun = systime(1) ; pour key_performance 83 @common 84 ;; 85 ;--------------------------------------------------------- 86 @cm_4mesh 87 IF NOT keyword_set(key_forgetold) THEN BEGIN 88 @updatenew 89 ENDIF 85 90 ;------------------------------------------------------------ 86 91 ; le masque est donne ou il faut prendre tmask? … … 92 97 ny = taille[2] 93 98 ; 94 ;------------------------------------------------------------ 95 if n_elements(periodique) EQ 0 then periodique = keyword_set(key_periodique) 96 if keyword_set(key_periodi que) and keyword_set(periodique) $97 AND NOT keyword_set( regulier) then BEGIN99 IF n_elements(keep_cont) EQ 0 THEN keep_cont = 1-key_irregular 100 ;------------------------------------------------------------ 101 if keyword_set(key_periodic)*(nx EQ jpi) $ 102 AND NOT keyword_set(basic) then BEGIN 98 103 msk = [msk, msk[0, *]] 99 104 nx = nx+1 … … 138 143 print, 'temps triangule: trouver pts_downward', systime(1)-tempdeux 139 144 ; 140 if NOT keyword_set(regulier) then begin145 if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 141 146 tempdeux = systime(1) ; pour key_performance =2 142 147 ;2 points terre en diagonale montante avec 2 points mer sur la diagonale descendante … … 156 161 print, 'temps triangule: trouver coindesc', systime(1)-tempdeux 157 162 ; 158 endif163 ENDIF 159 164 ; 160 165 if n_elements(pts_downward) EQ 1 then BEGIN … … 179 184 pts_downward =different(pts_downward,derniere_colonne ) 180 185 pts_downward =different(pts_downward,derniere_ligne ) 181 if NOT keyword_set(regulier) then begin186 if (NOT keyword_set(basic)) OR keyword_set(coinmonte) OR keyword_set(coindescend) then begin 182 187 if coinmont[0] NE -1 then begin 183 188 coinmont =different(coinmont,derniere_colonne ) … … 208 213 ; ca devrait aller mieux dans les prochaines versions d''IDL... 209 214 ; 210 tempdeux = systime(1) ; pour key_performance =2 215 if (NOT keyword_set(basic)) AND (NOT keyword_set(keep_cont)) then begin 216 tempdeux = systime(1) ; pour key_performance =2 211 217 ; on enleve les rectangles qui sont entierement dans la terre 212 recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1)213 IF testvar(var = key_performance) EQ 2 THEN $214 print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux218 recdsterre = where((1-msk)*(1-shift(msk, -1, 0))*(1-shift(msk, 0, -1))*(1-shift(msk, -1, -1)) EQ 1) 219 IF testvar(var = key_performance) EQ 2 THEN $ 220 print, 'temps triangule: tous les recdsterre', systime(1)-tempdeux 215 221 216 222 ; en attendant une version qui marche parfaitement, on est contraint … … 219 225 ; commun. 220 226 ; t1 = systime(1) 221 indice = intarr(nx, ny)222 trimask = intarr(nx, ny)223 trimask[0:nx-2, 0:ny-2] = 1224 IF recdsterre[0] NE -1 then BEGIN225 tempdeux = systime(1); pour key_performance =2226 indice[recdsterre] = 1227 if NOT keyword_set(regulier) then begin227 indice = intarr(nx, ny) 228 trimask = intarr(nx, ny) 229 trimask[0:nx-2, 0:ny-2] = 1 230 IF recdsterre[0] NE -1 then BEGIN 231 tempdeux = systime(1) ; pour key_performance =2 232 indice[recdsterre] = 1 233 ; if NOT keyword_set(basic) then begin 228 234 vire1 = 0 229 235 vire2 = 0 … … 247 253 IF testvar(var = key_performance) EQ 2 THEN $ 248 254 print, 'temps triangule: trier les recdsterre', systime(1)-tempdeux 255 ; endif 256 indice[*, ny-1] = 1 ; la deriere colonne te la derniere ligne 257 indice[nx-1, *] = 1 ; ne peuvent definir de rectangle. 258 ; 259 tempdeux = systime(1) ; pour key_performance =2 260 recgarde = where(indice EQ 0) 261 ; on recupere les numeros des triangles que l'' on va garder 262 trigarde = 2*[recgarde-recgarde/nx] 263 trigarde = transpose(temporary(trigarde)) 264 trigarde = [trigarde, trigarde+1] 265 ; 266 triang = triang[*, temporary(trigarde[*])] 267 IF testvar(var = key_performance) EQ 2 THEN $ 268 print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux 249 269 endif 250 indice[*, ny-1] = 1 ; la deriere colonne te la derniere ligne251 indice[nx-1, *] = 1 ; ne peuvent definirde rectangle.252 ;253 tempdeux = systime(1) ; pour key_performance =2254 recgarde = where(indice EQ 0)255 ; on recupere les numeros des triangles que l'' on va garder256 trigarde = 2*[recgarde-recgarde/nx]257 trigarde = [trigarde, trigarde+1]258 trigarde = trigarde[sort(trigarde)]259 ;260 triang = triang[*, trigarde]261 IF testvar(var = key_performance) EQ 2 THEN $262 print, 'temps triangule: virer les triangle de la liste', systime(1)-tempdeux263 270 endif 264 271 ; print, 'temps tri triangles', systime(1)-t1 265 272 ;------------------------------------------------------------ 266 ; quand key_periodi queeq 1, triang est une liste d''indice d'un273 ; quand key_periodic eq 1, triang est une liste d''indice d'un 267 274 ; tableau qui a une colonne de trop. 268 275 ; il faut ramener ca a la matrice initiale en mettant les indivces de … … 270 277 ;------------------------------------------------------------ 271 278 tempdeux = systime(1) ; pour key_performance =2 272 if keyword_set(key_periodi que) and keyword_set(periodique) $273 AND NOT keyword_set( regulier) then BEGIN279 if keyword_set(key_periodic)*(nx-1 EQ jpi) $ 280 AND NOT keyword_set(basic) then BEGIN 274 281 indicey = triang/nx 275 282 indicex = triang-indicey*nx … … 302 309 303 310 ;------------------------------------------------------------ 304 if arg_present(coinmonte) THEN coinmonte = coinmont ELSE cointerremont= coinmont305 if arg_present(coindescend) THEN coindescend = coindesc ELSE cointerredesc= coindesc306 ; 307 ; 308 ;------------------------------------------------------------ 309 311 if keyword_set(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont 312 if keyword_set(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc 313 ;------------------------------------------------------------ 314 IF NOT keyword_set(key_forgetold) THEN BEGIN 315 @updateold 316 ENDIF 310 317 311 318 IF keyword_set(key_performance) THEN print, 'temps triangule', systime(1)-tempsun -
trunk/ToBeReviewed/TRIANGULATION/triangule_e.pro
r27 r29 32 32 ;------------------------------------------------------------ 33 33 FUNCTION triangule_e, maskentree, COINMONTE = coinmonte, COINDESCEND = coindescend $ 34 , SHIFTED = shifted, REGULIER = regulier, PERIODIQUE = periodique 34 , SHIFTED = shifted, BASIC = basic 35 ;--------------------------------------------------------- 36 @cm_4mesh 37 IF NOT keyword_set(key_forgetold) THEN BEGIN 38 @updatenew 39 ENDIF 40 ;--------------------------------------------------------- 35 41 tempsun = systime(1) ; pour key_performance 36 @common37 ;;38 42 ;------------------------------------------------------------ 39 43 ; le masque est donne ou il faut prendre tmask? … … 45 49 ny = sizem[2] 46 50 ;------------------------------------------------------------ 47 if n_elements(periodique) EQ 0 then periodique = keyword_set(key_periodique) 48 if keyword_set(key_periodique) and keyword_set(periodique) $ 49 AND NOT keyword_set(regulier) then BEGIN 51 if keyword_set(key_periodic)*(nx EQ jpi) $ 52 AND NOT keyword_set(basic) then BEGIN 50 53 msk = [msk, msk[0, *]] 51 54 nx = nx+1 … … 132 135 ; ; 133 136 ;------------------------------------------------------------ 134 ; quand key_periodi queeq 1, triang est une liste d''indice d'un137 ; quand key_periodic eq 1, triang est une liste d''indice d'un 135 138 ; tableau qui a une colonne de trop. 136 139 ; il faut ramener ca a la matrice initiale en mettant les indivces de … … 138 141 ;------------------------------------------------------------ 139 142 tempdeux = systime(1) ; pour key_performance =2 140 if keyword_set(key_periodi que) and keyword_set(periodique) $141 AND NOT keyword_set( regulier) then BEGIN143 if keyword_set(key_periodic)*(nx-1 EQ jpi) $ 144 AND NOT keyword_set(basic) then BEGIN 142 145 indicey = triang/nx 143 146 indicex = triang-indicey*nx … … 170 173 171 174 ;------------------------------------------------------------ 172 ; if arg_present(coinmonte) THEN coinmonte = coinmont ELSE cointerremont= coinmont173 ; if arg_present(coindescend) THEN coindescend = coindesc ELSE cointerredesc= coindesc175 ; if arg_present(coinmonte) THEN coinmonte = coinmont ELSE twin_corners_up = coinmont 176 ; if arg_present(coindescend) THEN coindescend = coindesc ELSE twin_corners_dn = coindesc 174 177 ; 175 ; 178 ; IF NOT keyword_set(key_forgetold) THEN BEGIN 179 ; @updateold 180 ; ENDIF 176 181 ;------------------------------------------------------------ 177 182
Note: See TracChangeset
for help on using the changeset viewer.