Changeset 125 for trunk/SRC


Ignore:
Timestamp:
07/06/06 16:10:25 (18 years ago)
Author:
pinsard
Message:

improvements of Interpolation/*.pro header

Location:
trunk/SRC/Interpolation
Files:
22 edited

Legend:

Unmodified
Added
Removed

  • trunk/SRC/Interpolation/angle.pro

    r121 r125  
    11;--------- 
    22;+ 
    3 ; @file_comments north stereographic polar projection 
     3; @file_comments  
     4; north stereographic polar projection 
    45; 
    56; @param plam {in}{required} 
    67; 
    7 ; @param pphi {in}{required}  
     8; @param pphi {in}{required} 
    89; 
    910; @keyword DOUBLE {default=0} use double precision (default is float) 
    1011; 
    1112; @returns 
    12 ;       gsinu,gcosu : sinus and cosinus of the angle  
    13 ;       gsinv,gcosv   between north-south direction  
     13;       gsinu,gcosu : sinus and cosinus of the angle 
     14;       gsinv,gcosv   between north-south direction 
    1415;       gsint,gcost   and the j-direction of the mesh 
    1516; 
    16 ; @restrictions to compute the lateral boundary conditions, we assume 
    17 ; that:  
     17; @restrictions 
     18; to compute the lateral boundary conditions, we assume that: 
    1819;     (1) the first line is similar to the second line 
    19 ;       =>    gcosu[*, 0] = gcosu[*, 1]  
    20 ;       =>    gsinu[*, 0] = gsinu[*, 1]  
     20;       =>    gcosu[*, 0] = gcosu[*, 1] 
     21;       =>    gsinu[*, 0] = gsinu[*, 1] 
    2122;     (2) the grid follows OPA x periodicity rule, first column is 
    2223;     equal to the next to last column 
    … …  
    2627; 
    2728; @history 
    28 ;   -------------- 
    2929;       Original :  96-07 (O. Marti) 
    3030;                   98-06 (G. Madec) 
    31 ;       Feb 2005: IDL adaptation S. Masson  
     31;       Feb 2005: IDL adaptation S. Masson 
    3232; 
    3333; @version $Id$ 
    3434; 
    3535;- 
    36 ;--------- 
    3736; 
    3837FUNCTION fsnspp, plam, pphi, DOUBLE = double 
    … …  
    4746    a = 2. * tan( !pi/4. - !pi/180.*float(pphi)/2. ) 
    4847    x = cos( !pi/180.*float(plam) ) * a 
    49     y = sin( !pi/180.*float(plam) ) * a     
     48    y = sin( !pi/180.*float(plam) ) * a 
    5049  ENDELSE 
    5150  RETURN, {x:x, y:y} 
    … …  
    150149  znpv01 = -1; free memory 
    151150  IF keyword_set(double) THEN zmnpvt = 1.e-14 > zmnpvt $ 
    152   ELSE zmnpvt = 1.e-6 > zmnpvt  
     151  ELSE zmnpvt = 1.e-6 > zmnpvt 
    153152;   ... j-direction: f-point segment direction (u-point) 
    154153  zxffu = znpf00.x - znpf01.x 
    … …  
    157156  znpf01 = -1; free memory 
    158157  IF keyword_set(double) THEN zmnpfu = 1.e-14 > zmnpfu $ 
    159   ELSE zmnpfu = 1.e-6 > zmnpfu  
     158  ELSE zmnpfu = 1.e-6 > zmnpfu 
    160159;   ... i-direction: f-point segment direction (v-point) 
    161160  zxffv = znpf00.x - znpf10.x 
    162   zyffv = znpf00.y - znpf10.y  
     161  zyffv = znpf00.y - znpf10.y 
    163162  znpf00 = -1 &  znpf10 = -1; free memory 
    164163  zmnpfv = sqrt ( temporary(znnpv) * ( zxffv*zxffv + zyffv*zyffv )  ) 
    165164  IF keyword_set(double) THEN zmnpfv = 1.e-14 > zmnpfv $ 
    166   ELSE zmnpfv = 1.e-6 > zmnpfv  
     165  ELSE zmnpfv = 1.e-6 > zmnpfv 
    167166;   ... cosinus and sinus using scalar and vectorial products 
    168167  gsint = ( znpt.x*zyvvt - znpt.y*zxvvt ) / zmnpvt 
    … …  
    193192; ================================ 
    194193; 
    195   gcost[*, 0] = gcost[*, 1]  
    196   gsint[*, 0] = gsint[*, 1]  
    197   gcosu[*, 0] = gcosu[*, 1]  
    198   gsinu[*, 0] = gsinu[*, 1]  
     194  gcost[*, 0] = gcost[*, 1] 
     195  gsint[*, 0] = gsint[*, 1] 
     196  gcosu[*, 0] = gcosu[*, 1] 
     197  gsinu[*, 0] = gsinu[*, 1] 
    199198  gcosv[0, *] = gcosv[jpj-2, *] 
    200199  gsinv[0, *] = gsinv[jpj-2, *] 

  • trunk/SRC/Interpolation/clickincell.pro

    r121 r125  
    11;+ 
    2 ; @file_comments click on a map and find in which cell the click was 
     2; @file_comments  
     3; click on a map and find in which cell the click was 
    34; 
    45; @categories finding where is a point on a grid 
    … …  
    910;     points). 
    1011; 
    11 ; @keyword /DRAWCELL to draw the cell in which we clicked 
     12; @keyword DRAWCELL to draw the cell in which we clicked 
    1213; 
    1314; @keyword COLOR  the color used to draw the cells (Clicking one more 
    1415;     time in the same cell will draw the cell with the white color) 
    1516; 
    16 ; @keyword /ORIGINAL to get the position of the cell regarding the original 
     17; @keyword ORIGINAL to get the position of the cell regarding the original 
    1718;     grid (with no key_shift, ixminmesh, iyminmesh...) 
    1819; 
    19 ; @keyword /IJ see outpus 
     20; @keyword IJ see outpus 
    2021; 
    2122; @keyword _EXTRA to pass extra keywords to inquad and plot (when /drawcell) 
    … …  
    2526;     is in memory in the variable of the common. If /ij keyword is 
    2627;     activated give 2D array (2, n) which are the i,j position of the 
    27 ;     n selected cells.  
     28;     n selected cells. 
    2829; 
    2930; @uses common.pro 
    3031; 
    31 ; @examples  
     32; @examples 
    3233; 
    33 ;   IDL> res = clickincell() 
     34; IDL> res = clickincell() 
    3435;     Click with the left button to select a cell. Clicking one more 
    3536;     time in the same cell remove the cell from the selection. 
    36 ;     Click on the right button to quit.   
     37;     Click on the right button to quit. 
    3738; 
    38 ;   IDL> plt, findgen(jpi,jpj),/nodata,map=[90,0,0],/ortho 
    39 ;   IDL> print, clickincell(/draw,color=150,/xy) 
     39; IDL> plt, findgen(jpi,jpj),/nodata,map=[90,0,0],/ortho 
     40; IDL> print, clickincell(/draw,color=150,/xy) 
    4041; 
    4142; @history 
    … …  
    4445; 
    4546; @version $Id$ 
    46 ; 
    4747; 
    4848;- 
    … …  
    117117          cellnum = [cellnum, cell] 
    118118          selected = [selected, 1] 
    119           already = n_elements(selected)-1  
     119          already = n_elements(selected)-1 
    120120        ENDIF ELSE selected[already] = 1-selected[already] 
    121121        IF keyword_set(drawcell) THEN BEGIN 
    … …  
    129129      2:                        ; middle button 
    130130      ELSE: 
    131     ENDCASE     
     131    ENDCASE 
    132132; get mousse position on the reference map 
    133133outwhile: 
    … …  
    136136; 
    137137  good = where(selected NE 0) 
    138   IF good[0] EQ -1 THEN RETURN, -1  
     138  IF good[0] EQ -1 THEN RETURN, -1 
    139139; 
    140140  cellnum = cellnum[good] 
    … …  
    183183  ENDIF 
    184184; 
    185   ncell = n_elements(xx)  
     185  ncell = n_elements(xx) 
    186186  IF keyword_set(ij) THEN $ 
    187187    RETURN, [reform(xx, 1, ncell, /over) $ 
    188              , reform(yy, 1, ncell, /over)]  
     188             , reform(yy, 1, ncell, /over)] 
    189189; 
    190190  IF keyword_set(original) THEN RETURN, xx+jpiglo*yy $ 
    191   ELSE RETURN, xx+jpi*yy  
    192 END  
     191  ELSE RETURN, xx+jpi*yy 
     192END 

  • trunk/SRC/Interpolation/compute_fromirr_bilinear_weigaddr.pro

    r121 r125  
    11;+ 
    2 ; @file_comments compute the weight and address needed to interpolate data from 
    3 ;          an "irregular 2D grid" (defined as a grid made of quadrilateral cells) 
    4 ;          to any grid using the bilinear method 
    5 ;    
     2; @file_comments  
     3; compute the weight and address needed to interpolate data from 
     4; an "irregular 2D grid" (defined as a grid made of quadrilateral cells) 
     5; to any grid using the bilinear method 
     6; 
    67; @categories interpolation 
    78; 
    8 ; @param olonin {in}{required} longitudeof the input data  
    9 ; @param olat   {in}{required} latitude of the input data  
    10 ; @param omsk   {in}{required} land/se mask of the input data  
    11 ; @param alonin {in}{required} longitude of the output data  
    12 ; @param alat   {in}{required} latitude of the output data  
    13 ; @param amsk   {in}{required} land/sea mask of the output data  
     9; @param olonin {in}{required} longitudeof the input data 
     10; @param olat   {in}{required} latitude of the input data 
     11; @param omsk   {in}{required} land/se mask of the input data 
     12; @param alonin {in}{required} longitude of the output data 
     13; @param alat   {in}{required} latitude of the output data 
     14; @param amsk   {in}{required} land/sea mask of the output data 
    1415; 
    1516; @param weig {out} 
    … …  
    2122; 
    2223; @restrictions 
    23 ;  -  the input grid must be an "irregular 2D grid", defined as a grid made  
    24 ;     of quadrilateral cells which corners positions are defined with olonin and olat 
     24;  -  the input grid must be an "irregular 2D grid", defined as a grid made 
     25;  of quadrilateral cells which corners positions are defined with olonin and olat 
    2526;  -  We supposed the data are located on a sphere, with a periodicity along 
    26 ;     the longitude 
     27;  the longitude 
    2728;  -  to perform the bilinear interpolation within quadrilateral cells, we 
    28 ;     first morph the cell into a square cell and then compute the bilinear 
    29 ;     interpolation. 
    30 ;  -  if some corners of the cell are land points, their weight is set to 0  
    31 ;     and the weight is redistributed on the remaining "water" corners 
     29;  first morph the cell into a square cell and then compute the bilinear 
     30;  interpolation. 
     31;  -  if some corners of the cell are land points, their weight is set to 0 
     32;  and the weight is redistributed on the remaining "water" corners 
    3233;  -  points located out of the southern and northern boundaries or in cells 
    33 ;    containing only land points are set the the same value as their closest neighbor l;  
     34;  containing only land points are set the the same value as their closest neighbor l 
     35; 
    3436; @history 
    35 ;  June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr)  
    36 ;  
     37;  June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    3738; 
    3839; @version $Id$ 
    … …  
    4041;- 
    4142; 
    42 ;---------------------------------------------------------- 
    43 ;---------------------------------------------------------- 
    44 ; 
    4543PRO compute_fromirr_bilinear_weigaddr, olonin, olat, omsk, alonin, alat, amsk, weig, addr 
    4644; 
    47   compile_opt idl2, strictarrsubs  
     45  compile_opt idl2, strictarrsubs 
    4846; 
    4947  jpia = (size(alonin, /dimensions))[0] 
    … …  
    5553; longitude, between 0 and 360 
    5654  alon = alonin MOD 360 
    57   out = where(alon LT 0)  
     55  out = where(alon LT 0) 
    5856  IF out[0] NE -1 THEN alon[out] = alon[out]+360 
    5957  olon = olonin MOD 360 
    60   out = where(olon LT 0)  
     58  out = where(olon LT 0) 
    6159  IF out[0] NE -1 THEN olon[out] = olon[out]+360 
    6260; 
    6361; we work only on the water points 
    6462  owater = where(omsk EQ 1) 
    65   nowater = n_elements(owater)  
     63  nowater = n_elements(owater) 
    6664  awater = where(amsk EQ 1) 
    67   nawater = n_elements(awater)  
     65  nawater = n_elements(awater) 
    6866; 
    6967; define all cells that have corners located at olon, olat 
    … …  
    8785; list the cells that have at least 1 ocean point as corner 
    8886  good = where(total(omsk[celladdr], 1) GT 0) 
    89 ; keep only those cells  
     87; keep only those cells 
    9088  celladdr = celladdr[*, temporary(good)] 
    9189; 
    … …  
    110108    yy = alat[awater[n]] 
    111109; 1) we reduce the number of ocean cells for which we will try to know if 
    112 ; xx,yy is inside.  
     110; xx,yy is inside. 
    113111    CASE 1 OF 
    114112; if we are near the norh pole 
    … …  
    119117      END 
    120118; if we are near the longitude periodicity area 
    121       xx LE delta OR xx GE (360-delta):BEGIN  
     119      xx LE delta OR xx GE (360-delta):BEGIN 
    122120        lat1 = yy-delta 
    123121        lat2 = yy+delta 
    … …  
    157155      IF ind NE -1 THEN BEGIN 
    158156        ind = good[ind] 
    159 ; we keep its address  
     157; we keep its address 
    160158        foraddr[n] = ind 
    161159; now, we morph the quadrilateral ocean cell into the reference square (0 -> 1) 
    … …  
    173171                                    , xcell[2, ind], ycell[2, ind] $ 
    174172                                    , xcell[3, ind], ycell[3, ind], xx, yy) 
    175         ENDELSE  
     173        ENDELSE 
    176174; take care of rounding errors... 
    177175        zero = where(abs(xy) LT 1e-4) 

  • trunk/SRC/Interpolation/compute_fromreg_bilinear_weigaddr.pro

    r121 r125  
    11;+ 
    2 ; @file_comments compute the weight and address needed to interpolate data from a 
    3 ;          "regular grid" to any grid using the bilinear method 
    4 ;    
     2; @file_comments  
     3; compute the weight and address needed to interpolate data from a 
     4; "regular grid" to any grid using the bilinear method 
     5; 
    56; @categories interpolation 
    67; 
    7 ; @param alonin {in}{required} longitudeof the input data  
    8 ; @param alatin  {in}{required} latitude of the input data  
    9 ; @param olonin {in}{required} longitude of the output data  
    10 ; @param olat  {in}{required} latitude of the output data  
    11 ; 
    12 ; @keyword     /NONORTHERNLINE activate if you don't whant to take into 
     8; @param alonin {in}{required} longitudeof the input data 
     9; @param alatin  {in}{required} latitude of the input data 
     10; @param olonin {in}{required} longitude of the output data 
     11; @param olat  {in}{required} latitude of the output data 
     12; 
     13; @keyword NONORTHERNLINE activate if you don't want to take into 
    1314;          account the northen line of the input data when perfoming the 
    14 ; @keyword     /NOSOUTHERNLINE activate if you don't whant to take into 
     15; @keyword NOSOUTHERNLINE activate if you don't want to take into 
    1516;          account the southern line of the input data when perfoming the 
    1617;          interpolation. 
    … …  
    2425; 
    2526; @restrictions 
    26 ;  -  the input grid must be a "regular grid", defined as a grid for which each 
    27 ;     lontitudes lines have the same latitude and each latitudes columns have the 
    28 ;     same longitude. 
    29 ;  -  We supposed the data are located on a sphere, with a periodicity along 
    30 ;     the longitude. 
    31 ;  -  points located out of the southern and northern boundaries are interpolated 
    32 ;     using a linear interpolation only along the longitudinal direction. 
    33 ;  
     27;  - the input grid must be a "regular grid", defined as a grid for which each 
     28;  lontitudes lines have the same latitude and each latitudes columns have the 
     29;  same longitude. 
     30;  - We supposed the data are located on a sphere, with a periodicity along 
     31;  the longitude. 
     32;  - points located out of the southern and northern boundaries are interpolated 
     33;  using a linear interpolation only along the longitudinal direction. 
     34; 
    3435; @history 
    35 ;  November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr)  
    36 ;  
     36;  November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr) 
     37; 
    3738; @version $Id$ 
    3839; 
    3940;- 
    40 ; 
    41 ;---------------------------------------------------------- 
    42 ;---------------------------------------------------------- 
    4341; 
    4442PRO compute_fromreg_bilinear_weigaddr, alonin, alatin, olonin, olat, weig, addr $ 
    … …  
    6159  IF maxalon-minalon GE 360. THEN stop 
    6260; alon must be monotonically increasing 
    63   IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN BEGIN  
     61  IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN BEGIN 
    6462    shiftx = -(where(alon EQ min(alon)))[0] 
    6563    alon = shift(alon, shiftx) 
    6664    IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN stop 
    6765  ENDIF ELSE shiftx = 0 
    68 ; for longitude periodic bondary condition we add the fist 
    69 ; column on the right side of the array and  
     66; for longitude periodic boundary condition we add the fist 
     67; column on the right side of the array and 
    7068  alon = [alon, alon[0]+360.] 
    7169  jpia = jpia+1L 
    … …  
    7674  IF array_equal(sort(alat), lindgen(jpja)) NE 1 THEN stop 
    7775; 
    78   if keyword_set(nonorthernline) then BEGIN  
     76  if keyword_set(nonorthernline) then BEGIN 
    7977    jpja = jpja - 1L 
    8078    alat = alat[0: jpja-1L] 
    8179  ENDIF 
    82   if keyword_set(nosouthernline) then BEGIN  
     80  if keyword_set(nosouthernline) then BEGIN 
    8381    alat = alat[1: jpja-1L] 
    8482    jpja = jpja - 1L 
    … …  
    103101; if the ocean point is out of the atm grid, we use closest neighbor 
    104102; interpolation 
    105 ;  
     103; 
    106104; for each T point of oce grid, we find in which armospheric cell it is 
    107105; located. 
    … …  
    113111; for longitude, each ocean points must be located in atm cell. 
    114112  IF (where(pos[0, *] EQ -1))[0] NE -1 THEN stop 
    115 ; no ocean point should be located westward of the left bondary of the 
    116 ; atm cell in which it is supposed to be located  
     113; no ocean point should be located westward of the left boundary of the 
     114; atm cell in which it is supposed to be located 
    117115  IF total(olon LT alon[pos[0, *]]) NE 0 THEN stop 
    118 ; no ocean point should be located eastward of the right bondary of the 
    119 ; atm cell in which it is supposed to be located  
     116; no ocean point should be located eastward of the right boundary of the 
     117; atm cell in which it is supposed to be located 
    120118  IF total(olon GT alon[pos[0, *]+1]) NE 0 THEN stop 
    121119; 
    … …  
    124122; we change the coordinates of each ocean points to fit into a 
    125123; rectangle defined by: 
    126 ;  
     124; 
    127125;  y2 *------------* 
    128126;     |            | 
    … …  
    143141  IF max(indx) GT jpia-2 THEN stop ; checks... 
    144142  IF max(indy) GT jpja-2 THEN stop ; checks... 
    145 ; x coordinates of the atm cell  
     143; x coordinates of the atm cell 
    146144  x1 = alon[indx] 
    147145  x2 = alon[indx+1] 
    … …  
    149147  divi = temporary(x2)-x1 
    150148  glamnew = (olon-x1)/temporary(divi) 
    151   x1 = -1 ; free memory  
    152   olon = -1 ; free memory  
    153 ; y coordinates of the atm cell  
     149  x1 = -1 ; free memory 
     150  olon = -1 ; free memory 
     151; y coordinates of the atm cell 
    154152  y1 = alat[indy] 
    155153  y2 = alat[indy+1]             ; 
    … …  
    159157  IF zero[0] NE -1 THEN divi[zero] = 1. 
    160158  gphinew = (olat-y1)/temporary(divi) 
    161   y1 = -1 ; free memory  
     159  y1 = -1 ; free memory 
    162160; checks... 
    163161  IF min(glamnew) LT 0 THEN stop 
    … …  
    166164; weight and address array used for bilinear interpolation. 
    167165  xaddr = lonarr(4, jpio*jpjo) 
    168   xaddr[0, *] = indx   
     166  xaddr[0, *] = indx 
    169167  xaddr[1, *] = indx + 1L 
    170168  xaddr[2, *] = indx + 1L 
    171   xaddr[3, *] = indx   
    172 ;  
     169  xaddr[3, *] = indx 
     170; 
    173171  yaddr = lonarr(4, jpio*jpjo) 
    174172  yaddr[0, *] = indy 
    … …  
    181179  weig[1, *] =     glamnew  * (1.-gphinew) 
    182180  weig[2, *] =     glamnew  *     gphinew 
    183   weig[3, *] = (1.-glamnew) *     gphinew  
     181  weig[3, *] = (1.-glamnew) *     gphinew 
    184182; free memory 
    185183  gphinew = -1 
    … …  
    190188    ybad = olat[bad] 
    191189; the ocean points that are not located into an atm cell should be 
    192 ; located northward of the northern boundary of the atm grid  
    193 ;      or southward of the southern boundary of the atm grid  
     190; located northward of the northern boundary of the atm grid 
     191;      or southward of the southern boundary of the atm grid 
    194192    IF total(ybad GE min(alat) AND ybad LE max(alat)) GE 1 THEN stop 
    195193; 
    … …  
    223221  IF revy EQ 1 THEN yaddr = jpja - 1L - temporary(yaddr) 
    224222;                         ; 
    225   addr = temporary(yaddr)*jpia + temporary(xaddr)  
     223  addr = temporary(yaddr)*jpia + temporary(xaddr) 
    226224; 
    227225  return 

  • trunk/SRC/Interpolation/compute_fromreg_imoms3_weigaddr.pro

    r118 r125  
    11;+ 
    22; 
    3 ; @file_comments compute the weight and address neede to interpolate data from a 
    4 ;          "regular grid" to any grid using the imoms3 method 
    5 ;    
     3; @file_comments  
     4; compute the weight and address neede to interpolate data from a 
     5; "regular grid" to any grid using the imoms3 method 
     6; 
    67; @categories interpolation 
    78; 
    8 ; @param alonin {in}{required} longitude of the input data  
    9 ; @param alatin  {in}{required} latitude of the input data  
    10 ; @param olonin {in}{required} longitude of the output data  
    11 ; @param olat {in}{required} latitude of the output data  
    12 ; 
    13 ; @keyword /NONORTHERNLINE  
    14 ; @keyword /NOSOUTHERNLINE  
    15 ; activate if you don't whant to take into account the northen/southern line  
     9; @param alonin {in}{required} longitude of the input data 
     10; @param alatin {in}{required} latitude of the input data 
     11; @param olonin {in}{required} longitude of the output data 
     12; @param olat {in}{required} latitude of the output data 
     13; 
     14; @keyword NONORTHERNLINE 
     15; @keyword NOSOUTHERNLINE 
     16; activate if you don't want to take into account the northen/southern line 
    1617; of the input data when perfoming the interpolation. 
    1718; 
    … …  
    2425; 
    2526; @restrictions 
    26 ;  -  the input grid must be a "regular/rectangular grid", defined as a grid for 
    27 ;     which each lontitudes lines have the same latitude and each latitudes columns 
    28 ;     have the same longitude. 
     27;  - the input grid must be a "regular/rectangular grid", defined as a grid for 
     28;  which each lontitudes lines have the same latitude and each latitudes columns 
     29;  have the same longitude. 
    2930;  -  We supposed the data are located on a sphere, with a periodicity along 
    30 ;     the longitude. 
     31;  the longitude. 
    3132;  -  points located between the first/last 2 lines are interpolated 
    32 ;     using a imoms3 interpolation along the longitudinal direction and linear 
    33 ;     interpolation along the latitudinal direction 
    34 ;  -  points located out of the southern and northern boundaries are interpolated 
    35 ;     using a imoms3 interpolation only along the longitudinal direction. 
    36 ;  
     33;  using a imoms3 interpolation along the longitudinal direction and linear 
     34;  interpolation along the latitudinal direction 
     35;  - points located out of the southern and northern boundaries are interpolated 
     36;  using a imoms3 interpolation only along the longitudinal direction. 
     37; 
    3738; @history 
    38 ;  November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr)  
     39;  November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    3940;  March 2006: works for rectangular grids 
    4041; 
    … …  
    4950                                     , NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline 
    5051; 
    51   compile_opt idl2, strictarrsubs  
     52  compile_opt idl2, strictarrsubs 
    5253; 
    5354  alon = alonin 
    … …  
    6566  IF maxalon-minalon GE 360. THEN stop 
    6667; alon must be monotonically increasing 
    67   IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN BEGIN  
     68  IF array_equal(sort(alon), lindgen(jpia)) NE 1 THEN BEGIN 
    6869    shiftx = -(where(alon EQ min(alon)))[0] 
    6970    alon = shift(alon, shiftx) 
    … …  
    8990  IF total((step-step[0]) GE 1.e-6) NE 0 THEN noregy = 1 
    9091; 
    91   if keyword_set(nonorthernline) then BEGIN  
     92  if keyword_set(nonorthernline) then BEGIN 
    9293    jpja = jpja - 1L 
    9394    alat = alat[0: jpja-1L] 
    9495  ENDIF 
    95   if keyword_set(nosouthernline) then BEGIN  
     96  if keyword_set(nosouthernline) then BEGIN 
    9697    alat = alat[1: jpja-1L] 
    9798    jpja = jpja - 1L 
    … …  
    123124  indexlat = value_locate(alat, olat) 
    124125; 
    125 ; for the ocean points located below the atm line  
     126; for the ocean points located below the atm line 
    126127; jpja-2 and above the line 1 
    127128; for those points we can always find 16 neighbors 
    … …  
    131132  ilon = indexlon[short] 
    132133  ilat = indexlat[short] 
    133 ;  
    134   IF NOT keyword_set(noregy) THEN BEGIN  
     134; 
     135  IF NOT keyword_set(noregy) THEN BEGIN 
    135136    delta = alat[ilat+1L]-alat[ilat] 
    136137    IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop 
    … …  
    147148    d2 = (alat[ilat+1L]-olat[short])/delta 
    148149    IF min(d2, max = ma) LE 0 THEN stop 
    149     IF ma GT 1 THEN stop   
     150    IF ma GT 1 THEN stop 
    150151    wy2 = imoms3(temporary(d2)) 
    151152    d3 = (alat[ilat+2L]-olat[short])/delta 
    152153    IF min(d3, max = ma) LE 1 THEN stop 
    153     IF ma GT 2 THEN stop   
     154    IF ma GT 2 THEN stop 
    154155    wy3 = imoms3(temporary(d3)) 
    155   ENDIF ELSE BEGIN  
     156  ENDIF ELSE BEGIN 
    156157    nele = n_elements(short) 
    157158    wy0 = fltarr(nele) 
    … …  
    169170      wy3[i] = imoms3(2-newlat) 
    170171    ENDFOR 
    171   ENDELSE  
     172  ENDELSE 
    172173; 
    173174  mi = min(wy0+wy1+wy2+wy3, max = ma) 
    … …  
    175176  IF abs(ma-1) GE 1.e-6 THEN stop 
    176177; 
    177   IF NOT keyword_set(noregx) THEN BEGIN  
     178  IF NOT keyword_set(noregx) THEN BEGIN 
    178179    delta = alon[ilon]-alon[ilon-1L] 
    179180    IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop 
    … …  
    194195    d3 = (alon[ilon+2L]-olon[short])/delta 
    195196    IF min(d3, max = ma) LE 1 THEN stop 
    196     IF ma GT 2 THEN stop   
     197    IF ma GT 2 THEN stop 
    197198    wx3 = imoms3(temporary(d3)) 
    198   ENDIF ELSE BEGIN  
     199  ENDIF ELSE BEGIN 
    199200    nele = n_elements(short) 
    200201    wx0 = fltarr(nele) 
    … …  
    212213      wx3[i] = imoms3(2-newlon) 
    213214    ENDFOR 
    214   ENDELSE  
     215  ENDELSE 
    215216; 
    216217  mi = min(wx0+wx1+wx2+wx3, max = ma) 
    … …  
    220221; line 0 
    221222  xaddr[0, short] = ilon - 1L 
    222   xaddr[1, short] = ilon   
     223  xaddr[1, short] = ilon 
    223224  xaddr[2, short] = ilon + 1L 
    224225  xaddr[3, short] = ilon + 2L 
    … …  
    233234; line 1 
    234235  xaddr[4, short] = ilon - 1L 
    235   xaddr[5, short] = ilon   
     236  xaddr[5, short] = ilon 
    236237  xaddr[6, short] = ilon + 1L 
    237238  xaddr[7, short] = ilon + 2L 
    238   yaddr[4, short] = ilat   
    239   yaddr[5, short] = ilat   
    240   yaddr[6, short] = ilat   
    241   yaddr[7, short] = ilat   
     239  yaddr[4, short] = ilat 
     240  yaddr[5, short] = ilat 
     241  yaddr[6, short] = ilat 
     242  yaddr[7, short] = ilat 
    242243  weig[4, short] = wx0 * wy1 
    243244  weig[5, short] = wx1 * wy1 
    … …  
    246247; line 2 
    247248  xaddr[8, short] = ilon - 1L 
    248   xaddr[9, short] = ilon   
     249  xaddr[9, short] = ilon 
    249250  xaddr[10, short] = ilon + 1L 
    250251  xaddr[11, short] = ilon + 2L 
    … …  
    259260; line 3 
    260261  xaddr[12, short] = ilon - 1L 
    261   xaddr[13, short] = ilon   
     262  xaddr[13, short] = ilon 
    262263  xaddr[14, short] = ilon + 1L 
    263264  xaddr[15, short] = ilon + 2L 
    … …  
    275276  IF abs(ma-1) GE 1.e-6 THEN stop 
    276277; 
    277 ; for the ocean points located between the atm lines  
     278; for the ocean points located between the atm lines 
    278279; jpja-2 and jpja-1 or between the atm lines 0 and 1 
    279280; linear interpolation between line 1 and line 2 
    … …  
    285286; 
    286287    delta = alat[ilat+1L]-alat[ilat] 
    287     IF NOT keyword_set(noregy) THEN BEGIN  
     288    IF NOT keyword_set(noregy) THEN BEGIN 
    288289      IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop 
    289       delta = delta[0]  
    290     ENDIF  
     290      delta = delta[0] 
     291    ENDIF 
    291292; 
    292293    d1 = (alat[ilat   ]-olat[short])/delta 
    … …  
    296297    d2 = (alat[ilat+1L]-olat[short])/delta 
    297298    IF min(d2, max = ma) LE 0 THEN stop 
    298     IF ma GT 1 THEN stop   
     299    IF ma GT 1 THEN stop 
    299300    wy2 = 1.- temporary(d2) 
    300301; 
    … …  
    303304    IF abs(ma-1) GE 1.e-6 THEN stop 
    304305; but imoms3 along the longitude 
    305     IF NOT keyword_set(noregx) THEN BEGIN  
     306    IF NOT keyword_set(noregx) THEN BEGIN 
    306307      delta = alon[ilon]-alon[ilon-1L] 
    307308      IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop 
    … …  
    322323      d3 = (alon[ilon+2L]-olon[short])/delta 
    323324      IF min(d3, max = ma) LE 1 THEN stop 
    324       IF ma GT 2 THEN stop   
     325      IF ma GT 2 THEN stop 
    325326      wx3 = imoms3(temporary(d3)) 
    326     ENDIF ELSE BEGIN  
     327    ENDIF ELSE BEGIN 
    327328      nele = n_elements(short) 
    328329      wx0 = fltarr(nele) 
    … …  
    340341        wx3[i] = imoms3(2-newlon) 
    341342      ENDFOR 
    342     ENDELSE  
     343    ENDELSE 
    343344; 
    344345    mi = min(wx0+wx1+wx2+wx3, max = ma) 
    … …  
    347348; line 1 
    348349    xaddr[0, short] = ilon - 1L 
    349     xaddr[1, short] = ilon   
     350    xaddr[1, short] = ilon 
    350351    xaddr[2, short] = ilon + 1L 
    351352    xaddr[3, short] = ilon + 2L 
    352     yaddr[0, short] = ilat   
    353     yaddr[1, short] = ilat   
    354     yaddr[2, short] = ilat   
    355     yaddr[3, short] = ilat   
     353    yaddr[0, short] = ilat 
     354    yaddr[1, short] = ilat 
     355    yaddr[2, short] = ilat 
     356    yaddr[3, short] = ilat 
    356357    weig[0, short] = wx0 * wy1 
    357358    weig[1, short] = wx1 * wy1 
    … …  
    360361; line 2 
    361362    xaddr[4, short] = ilon - 1L 
    362     xaddr[5, short] = ilon   
     363    xaddr[5, short] = ilon 
    363364    xaddr[6, short] = ilon + 1L 
    364365    xaddr[7, short] = ilon + 2L 
    … …  
    379380; 
    380381; for the ocean points located below the line 0 
    381 ; Interpolation only along the longitude  
     382; Interpolation only along the longitude 
    382383; 
    383384  short = where(indexlat EQ -1) 
    … …  
    385386    ilon = indexlon[short] 
    386387; 
    387     IF NOT keyword_set(noregx) THEN BEGIN  
     388    IF NOT keyword_set(noregx) THEN BEGIN 
    388389      delta = alon[ilon]-alon[ilon-1L] 
    389390      IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop 
    … …  
    404405      d3 = (alon[ilon+2L]-olon[short])/delta 
    405406      IF min(d3, max = ma) LE 1 THEN stop 
    406       IF ma GT 2 THEN stop   
     407      IF ma GT 2 THEN stop 
    407408      wx3 = imoms3(temporary(d3)) 
    408     ENDIF ELSE BEGIN  
     409    ENDIF ELSE BEGIN 
    409410      nele = n_elements(short) 
    410411      wx0 = fltarr(nele) 
    … …  
    422423        wx3[i] = imoms3(2-newlon) 
    423424      ENDFOR 
    424     ENDELSE  
     425    ENDELSE 
    425426; 
    426427    mi = min(wx0+wx1+wx2+wx3, max = ma) 
    … …  
    429430; line 1 
    430431    xaddr[0, short] = ilon - 1L 
    431     xaddr[1, short] = ilon   
     432    xaddr[1, short] = ilon 
    432433    xaddr[2, short] = ilon + 1L 
    433434    xaddr[3, short] = ilon + 2L 
    434     yaddr[0:3, short] = 0  
     435    yaddr[0:3, short] = 0 
    435436    weig[0, short] = wx0 
    436437    weig[1, short] = wx1 
    … …  
    444445  ENDIF 
    445446; 
    446 ; for the ocean points located above jpia-1  
    447 ; Interpolation only along the longitude  
     447; for the ocean points located above jpia-1 
     448; Interpolation only along the longitude 
    448449; 
    449450  short = where(indexlat EQ jpja-1L) 
    … …  
    451452    ilon = indexlon[short] 
    452453; 
    453     IF NOT keyword_set(noregx) THEN BEGIN  
     454    IF NOT keyword_set(noregx) THEN BEGIN 
    454455      delta = alon[ilon]-alon[ilon-1L] 
    455456      IF max(abs(delta-delta[0])) GE 1.e-6 THEN stop 
    … …  
    470471      d3 = (alon[ilon+2L]-olon[short])/delta 
    471472      IF min(d3, max = ma) LE 1 THEN stop 
    472       IF ma GT 2 THEN stop   
     473      IF ma GT 2 THEN stop 
    473474      wx3 = imoms3(temporary(d3)) 
    474     ENDIF ELSE BEGIN  
     475    ENDIF ELSE BEGIN 
    475476      nele = n_elements(short) 
    476477      wx0 = fltarr(nele) 
    … …  
    488489        wx3[i] = imoms3(2-newlon) 
    489490      ENDFOR 
    490     ENDELSE  
     491    ENDELSE 
    491492; 
    492493    mi = min(wx0+wx1+wx2+wx3, max = ma) 
    … …  
    495496; line 1 
    496497    xaddr[0, short] = ilon-1L 
    497     xaddr[1, short] = ilon   
     498    xaddr[1, short] = ilon 
    498499    xaddr[2, short] = ilon+1L 
    499500    xaddr[3, short] = ilon+2L 
    … …  
    531532  IF revy EQ 1 THEN yaddr = jpja - 1L - temporary(yaddr) 
    532533;                        ; 
    533   addr = temporary(yaddr)*jpia+temporary(xaddr)  
     534  addr = temporary(yaddr)*jpia+temporary(xaddr) 
    534535; 
    535536  RETURN 

  • trunk/SRC/Interpolation/cutpar.pro

    r121 r125  
    11;+ 
    22; 
    3 ; @file_comments cut p parallelogram(s) into p*n^2 parallelograms 
     3; @file_comments  
     4; cut p parallelogram(s) into p*n^2 parallelograms 
    45; 
    56; @categories basic work 
    67; 
    78; @param x0 {in}{required} 
    8 ; @param y0 {in}{required}  
     9; @param y0 {in}{required} 
    910; @param x1 {in}{required} 
    10 ; @param y1 {in}{required}  
     11; @param y1 {in}{required} 
    1112; @param x2 {in}{required} 
    12 ; @param y2 {in}{required}  
     13; @param y2 {in}{required} 
    1314; @param x3 {in}{required} 
    14 ; @param y3 {in}{required}  
     15; @param y3 {in}{required} 
    1516; 1d arrays of p elements, giving the edge positions. The 
    16 ;       edges must be given as in plot to traw the parallelogram. (see 
     17;       edges must be given as in plot to draw the parallelogram. (see 
    1718;       example). 
    1819; 
    1920; @param n {in}{required} each parallelogram will be cutted in n^2 pieces 
    2021; 
    21 ; @keyword /endpoints see outputs 
     22; @keyword ENDPOINTS see outputs 
    2223; 
    23 ; @keyword /onsphere to specify that the points are located on a 
     24; @keyword ONSPHERE to specify that the points are located on a 
    2425;         sphere. In this case, x and y corresponds to longitude and 
    2526;         latitude in degrees. 
    2627; 
    2728; @returns 
    28 ;        - default: 3d array(2,n^2,p) giving the center position of each 
    29 ;        piece of the parallelograms 
    30 ;        - /endpoints: 3d array(2,(n+1)^2,p) giving the edge positions 
    31 ;        of each piece of the parallelograms 
     29;  - default: a 3d array(2,n^2,p) giving the center position of each 
     30;  piece of the parallelograms 
     31;  - if /ENDPOINTS : a 3d array(2,(n+1)^2,p) giving the edge positions 
     32;  of each piece of the parallelograms 
    3233; 
    3334; @uses cutsegment.pro 
    3435; 
    35 ; @examples  
     36; @examples 
    3637; 
    3738; IDL> x0 = [2,6,2] 
    … …  
    5657; 
    5758;- 
    58 FUNCTION cutpar, x0, y0, x1, y1, x2, y2, x3, y3, n, endpoints = endpoints, onsphere = onsphere 
     59FUNCTION cutpar, x0, y0, x1, y1, x2, y2, x3, y3, n, ENDPOINTS = endpoints, ONSPHERE = onsphere 
    5960; 
    6061  compile_opt idl2, strictarrsubs 
    … …  
    6667;   THEN stop; print, 'NOT a parallelogram' 
    6768; x0(npar) 
    68   npar = n_elements(x0)  
     69  npar = n_elements(x0) 
    6970; firstborder(2,n+keyword_set(endpoints),npar) 
    7071  firstborder = cutsegment(x0, y0, x1, y1, n $ 

  • trunk/SRC/Interpolation/cutsegment.pro

    r121 r125  
    11;+ 
    22; 
    3 ; @file_comments cut p segments into p*n equal parts 
     3; @file_comments  
     4; cut p segments into p*n equal parts 
    45; 
    56; @categories basic work 
    … …  
    1314; @param n {in}{required} the number of pieces we want to cut each segment 
    1415; 
    15 ; @keyword /endpoints see ouputs 
     16; @keyword ENDPOINTS see ouputs 
    1617; 
    17 ; @keyword /onsphere to specify that the points are located on a 
     18; @keyword ONSPHERE to specify that the points are located on a 
    1819;         sphere. In this case, x and y corresponds to longitude and 
    1920;         latitude in degrees. 
    2021; 
    2122; @returns 
    22 ;        default: a 3d array (2,n,p) that gives the coordinates of the 
    23 ;        middle of the cutted segments. 
    24 ;        if /endpoints, a 3d array (2,n+1,p) that gives the 
    25 ;        coordinates of the endpoints of the cutted segments. 
     23;  - default: a 3d array (2,n,p) that gives the coordinates of the 
     24;  middle of the cutted segments. 
     25;  - if /ENDPOINTS, a 3d array (2,n+1,p) that gives the 
     26;  coordinates of the endpoints of the cutted segments. 
    2627; 
    27 ; @examples  
     28; @examples 
    2829; 
    29 ;  IDL> x0=[2,5] 
    30 ;  IDL> y0=[5,1] 
    31 ;  IDL> x1=[9,3] 
    32 ;  IDL> y1=[1,8] 
    33 ;  IDL> res=cutsegment(x0, y0, x1, y1, 10) 
    34 ;  IDL> splot, [0,10], [0,10], xstyle = 1, ystyle = 1,/nodata 
    35 ;  IDL> oplot, [x0[0], x1[0]], [y0[0], y1[0]] 
    36 ;  IDL> oplot, [res[0,*,0]], [res[1,*,0]], color = 20, psym = 1, thick = 3 
    37 ;  IDL> oplot, [x0[1], x1[1]], [y0[1], y1[1]] 
    38 ;  IDL> oplot, [res[0,*,1]], [res[1,*,1]], color = 40, psym = 1, thick = 3 
     30; IDL> x0=[2,5] 
     31; IDL> y0=[5,1] 
     32; IDL> x1=[9,3] 
     33; IDL> y1=[1,8] 
     34; IDL> res=cutsegment(x0, y0, x1, y1, 10) 
     35; IDL> splot, [0,10], [0,10], xstyle = 1, ystyle = 1,/nodata 
     36; IDL> oplot, [x0[0], x1[0]], [y0[0], y1[0]] 
     37; IDL> oplot, [res[0,*,0]], [res[1,*,0]], color = 20, psym = 1, thick = 3 
     38; IDL> oplot, [x0[1], x1[1]], [y0[1], y1[1]] 
     39; IDL> oplot, [res[0,*,1]], [res[1,*,1]], color = 40, psym = 1, thick = 3 
    3940; 
    4041; @history 
    … …  
    4546; 
    4647;- 
    47 FUNCTION cutsegment, x0, y0, x1, y1, n, endpoints = endpoints, onsphere = onsphere 
     48FUNCTION cutsegment, x0, y0, x1, y1, n, ENDPOINTS = endpoints, ONSPHERE = onsphere 
    4849; 
    4950  compile_opt idl2, strictarrsubs 
    5051; 
    5152; number of segment 
    52   nseg = n_elements(x0)  
     53  nseg = n_elements(x0) 
    5354; number of point to find on each segment 
    54   n2find = n+keyword_set(endpoints)  
     55  n2find = n+keyword_set(endpoints) 
    5556; 
    5657  IF keyword_set(onsphere) THEN BEGIN 

  • trunk/SRC/Interpolation/extrapolate.pro

    r118 r125  
    11;+ 
    2 ; @file_comments extrapolate data (zinput) where maskinput eq 0 by filling  
     2; @file_comments  
     3; extrapolate data (zinput) where maskinput eq 0 by filling 
    34; step by step the coastline points with the mean value of the 8 neighbourgs. 
    45; 
    … …  
    2021FUNCTION extrapolate, zinput, maskinput, nb_iteration, x_periodic = x_periodic, MINVAL = minval, MAXVAL = maxval 
    2122; 
    22   compile_opt idl2, strictarrsubs  
     23  compile_opt idl2, strictarrsubs 
    2324; 
    2425; check the number of iteration used in the extrapolation. 
    … …  
    2829  ny = (size(zinput))[2] 
    2930; take care of the boundary conditions... 
    30 ;  
     31; 
    3132; for the x direction, we put 2 additional columns at the left and 
    32 ; right side of the array.  
     33; right side of the array. 
    3334; for the y direction, we put 2 additional lines at the bottom and 
    34 ; top side of the array.  
     35; top side of the array. 
    3536; These changes allow us to use shift function without taking care of 
    3637; the x and y periodicity. 
    … …  
    5455;--------------------------------------------------------------- 
    5556;--------------------------------------------------------------- 
    56 ; extrapolation  
     57; extrapolation 
    5758;--------------------------------------------------------------- 
    5859  sqrtinv = 1./sqrt(2) 
    5960; 
    6061  cnt = 1 
    61 ; When we look for the coast line, we don't whant to select the 
     62; When we look for the coast line, we don't want to select the 
    6263; borderlines of the array. -> we force the value of the mask for 
    6364; those lines. 
    … …  
    8687; we compute the weighted number of sea neighbourgs. 
    8788; those 4 neighbours have a weight of 1: 
    88 ;    *     
    89 ;   *+*         
    90 ;    *      
     89;    * 
     90;   *+* 
     91;    * 
    9192; 
    9293; those 4 neighbours have a weight of 1/sqrt(2): 
    9394; 
    9495;    * * 
    95 ;     +        
     96;     + 
    9697;    * * 
    9798; 
    … …  
    117118             +1./sqrt(2)*(z[nx2+1+coast]+z[nx2-1+coast] $ 
    118119                          +z[-nx2+1+coast]+z[-nx2-1+coast]) 
    119 ;     
     120; 
    120121    IF n_elements(minval) NE 0 THEN zcoast = minval > temporary(zcoast) 
    121122    IF n_elements(maxval) NE 0 THEN zcoast = temporary(zcoast) < maxval 
    … …  
    132133; 
    133134    cnt = cnt + 1 
    134 ; When we look for the coast line, we don't whant to select the 
     135; When we look for the coast line, we don't want to select the 
    135136; borderlines of the array. -> we force the value of the mask for 
    136137; those lines. 
    … …  
    146147; we return the original size of the array 
    147148;--------------------------------------------------------------- 
    148    
     149; 
    149150  return, z[1:nx, 1:ny] 
    150 END  
     151END 
    151152 

  • trunk/SRC/Interpolation/fromirr.pro

    r121 r125  
    11;+ 
    22; 
    3 ; @file_comments interpolate data from an irregular 2D grid to any 2D grid. 
    4 ;   Only 1 metod available = bilinear 
    5 ;    
     3; @file_comments  
     4; interpolate data from an irregular 2D grid to any 2D grid. 
     5;   Only 1 method available = bilinear 
     6; 
    67; @categories interpolation 
    78; 
    8 ;    @param method {in}{required} a string defining the interpolation method. must be 'bilinear' 
    9 ;    @param datain {in}{required} a 2D array the input data to interpolate 
    10 ;    @param lonin {in}{optional} a 2D array defining the longitude of the input data 
    11 ;    @param latin {in}{optional} a 2D array defining the latitude of the input data. 
    12 ;    @param mskin {in}{optional} a 2D array, the land-sea mask of the input data (1 on ocean, 0 on land) 
    13 ;    @param lonout {in}{optional} 1D or 2D array defining the longitude of the output data. 
    14 ;    @param latout {in}{optional} 1D or 2D array defining the latitude of the output data. 
    15 ;    @param mskout {in}{required} a 2D array, the land-sea mask of the ouput data (1 on ocean, 0 on land) 
     9; @param method {in}{required} a string defining the interpolation method. must be 'bilinear' 
     10; @param datain {in}{required} a 2D array the input data to interpolate 
     11; @param lonin {in}{optional} a 2D array defining the longitude of the input data 
     12; @param latin {in}{optional} a 2D array defining the latitude of the input data. 
     13; @param mskin {in}{optional} a 2D array, the land-sea mask of the input data (1 on ocean, 0 on land) 
     14; @param lonout {in}{optional} 1D or 2D array defining the longitude of the output data. 
     15; @param latout {in}{optional} 1D or 2D array defining the latitude of the output data. 
     16; @param mskout {in}{required} a 2D array, the land-sea mask of the ouput data (1 on ocean, 0 on land) 
    1617; 
    1718; @keyword WEIG (see ADDR) 
    … …  
    2728; @returns 2D array the interpolated data 
    2829; 
    29 ; @restrictions We supposed the data are located on a sphere, with a periodicity along 
    30 ;              the longitude. 
    31 ;              Note that the input data can contain the same cells several times 
    32 ;             (like ORCA grid near the north pole boundary) 
     30; @restrictions 
     31; We supposed the data are located on a sphere, with a periodicity along 
     32; the longitude. 
     33; Note that the input data can contain the same cells several times 
     34; (like ORCA grid near the north pole boundary) 
    3335; 
    3436; @examples 
    35 ;   
     37; 
    3638; IDL> tncep = fromirr('bilinear', topa, glamt, gphit, tmask[*,*,0], lonout, latout, mskout) 
    3739; 
    38 ;  or  
     40;  or 
    3941; 
    4042; IDL> t1ncep = fromirr('bilinear', topa, glamt, gphit, tmask[*,*,0], lonout, latout, mskout $ 
    … …  
    4446; 
    4547; @history 
    46 ;  June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr)  
    47 ;  
     48;  June 2006: Sebastien Masson (smasson\@lodyc.jussieu.fr) 
     49; 
    4850; @version $Id$ 
    4951; 
    … …  
    5557                  , WEIG = weig, ADDR = addr 
    5658; 
    57   compile_opt strictarr, strictarrsubs  
     59  compile_opt strictarr, strictarrsubs 
    5860; 
    5961;--------------- 
    … …  
    7577    CASE method OF 
    7678      'bilinear':compute_fromirr_bilinear_weigaddr, alon, alat, mskin, olon, olat, mskout, weig, addr 
    77       ELSE:BEGIN  
     79      ELSE:BEGIN 
    7880        print, ' unknown interpolation method... we stop' 
    7981        stop 

  • trunk/SRC/Interpolation/fromreg.pro

    r121 r125  
    11;+ 
    22; 
    3 ; @file_comments interpolate data from a "regular/rectangular grid" to any grid. 
    4 ;   2 metods availables: bilinear and imoms3  
    5 ;   A "regular/rectangular grid" is defined as a grid for which each lontitudes lines have  
     3; @file_comments  
     4; interpolate data from a "regular/rectangular grid" to any grid. 
     5;   2 methods availables: bilinear and imoms3 
     6;   A "regular/rectangular grid" is defined as a grid for which each lontitudes lines have 
    67;   the same latitude and each latitudes columns have the same longitude. 
    7 ;    
     8; 
    89; @categories interpolation 
    910; 
    10 ; @param method {in}{required}  a string defining the interpolation method.  
     11; @param method {in}{required}  a string defining the interpolation method. 
    1112;            must be 'bilinear' or 'imoms3' 
    1213; @param datain {in}{required}  a 2D array the input data to interpolate 
    … …  
    2627;     case, lonin, latin, lonout and latout are not necessary. 
    2728; 
    28 ; @keyword /NONORTHERNLINE  
    29 ; @keyword /NOSOUTHERNLINE  
    30 ; activate if you don't whant to take into account the northen/southern line  
     29; @keyword NONORTHERNLINE 
     30; @keyword NOSOUTHERNLINE 
     31; activate if you don't want to take into account the northen/southern line 
    3132; of the input data when perfoming the interpolation. 
    3233; 
    3334; @returns 2D array the interpolated data 
    3435; 
    35 ; @restrictions We supposed the data are located on a sphere, with a  
    36 ; periodicity along the longitude. 
     36; @restrictions 
     37; We supposed the data are located on a sphere, with a periodicity along the 
     38; longitude. 
    3739; 
    38 ; @examples   
    39 ;   
    40 ;  IDL> topa = fromreg('bilinear', tncep, xncep, yncep, glamt, gphit) 
     40; @examples 
    4141; 
    42 ;  or  
     42; IDL> topa = fromreg('bilinear', tncep, xncep, yncep, glamt, gphit) 
    4343; 
    44 ;  IDL> t1opa = fromreg('bilinear', t1ncep, xncep, yncep, glamt, gphit, WEIG = a, ADDR = b) 
    45 ;  IDL> help, a, b 
    46 ;  IDL> t2opa = fromreg('bilinear', t2ncep, xncep, WEIG = a, ADDR = b) 
     44;  or 
     45; 
     46; IDL> t1opa = fromreg('bilinear', t1ncep, xncep, yncep, glamt, gphit, WEIG = a, ADDR = b) 
     47; IDL> help, a, b 
     48; IDL> t2opa = fromreg('bilinear', t2ncep, xncep, WEIG = a, ADDR = b) 
    4749; 
    4850; @history 
    49 ;  November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr)  
     51;  November 2005: Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    5052; 
    5153; @version $Id$ 
    … …  
    6062                  , NOSOUTHERNLINE = nosouthernline 
    6163; 
    62   compile_opt idl2, strictarrsubs  
     64  compile_opt idl2, strictarrsubs 
    6365; 
    6466;--------------- 
    … …  
    8183      'bilinear':compute_fromreg_bilinear_weigaddr, alon, alat, olon, olat, weig, addr, NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline 
    8284      'imoms3':  compute_fromreg_imoms3_weigaddr,   alon, alat, olon, olat, weig, addr, NONORTHERNLINE = nonorthernline, NOSOUTHERNLINE = nosouthernline 
    83       ELSE:BEGIN  
     85      ELSE:BEGIN 
    8486        print, ' unknown interpolation method... we stop' 
    8587        stop 

  • trunk/SRC/Interpolation/get_gridparams.pro

    r121 r125  
    55;      and make sure it is 1D or 2D arrays 
    66; 
    7 ;   or  
     7;   or 
    88;   2) given longitude and latitude arrays get their dimensions and make 
    99;  sure they are 1D or 2D arrays 
    … …  
    1313; @examples 
    1414; 
    15 ; 1)  
     15; 1) 
    1616; IDL> get_gridparams, file, lonname, latname, lon, lat, jpi, jpj, n_dimensions 
    1717; 
    1818; or 
    1919; 
    20 ; 2)  
     20; 2) 
    2121; IDL> get_gridparams, lon, lat, jpi, jpj, n_dimensions 
    2222; 
    … …  
    4747;    and each latitudes should be the sae for all longitudes). 
    4848; 
    49 ; @keyword /DOUBLE use double precision to perform the computation 
     49; @keyword DOUBLE use double precision to perform the computation 
    5050; 
    5151; @examples 

  • trunk/SRC/Interpolation/imoms3.pro

    r118 r125  
    11;+ 
    2 ;  
    3 ;  
     2; 
     3; 
    44; @param xin {in}{required} 
    55; 
    … …  
    1313 
    1414  x = abs(xin) 
    15   y = fltarr(n_elements(x))  
     15  y = fltarr(n_elements(x)) 
    1616 
    1717  test1 = where(x LT 1) 

  • trunk/SRC/Interpolation/inquad.pro

    r121 r125  
    11;+ 
    2 ; @file_comments to find if an (x,y) point is in a quadrilateral (x1,x2,x3,x4) 
     2; @file_comments  
     3; to find if an (x,y) point is in a quadrilateral (x1,x2,x3,x4) 
    34; 
    45; @categories grid manipulation 
    … …  
    67; @param x {in}{required} 
    78; @param y {in}{required} 
    8 ;  the coordinates of the point we want to know where it is.  
    9 ;  Must be a scalar if /onsphere activated else can be scalar or array.  
     9;  the coordinates of the point we want to know where it is. 
     10;  Must be a scalar if /ONSPHERE activated else can be scalar or array. 
    1011; 
    1112; @param x1 {in}{required} 
    … …  
    1718; @param x4 {in}{required} 
    1819; @param y4 {in}{required} 
    19 ;  the coordinates of the quadrilateral given in the CLOCKWISE order.  
     20;  the coordinates of the quadrilateral given in the CLOCKWISE order. 
    2021;  Scalar or array. 
    2122; 
    22 ; @keyword /DOUBLE use double precision to perform the computation  
    23 ; 
    24 ; @keyword /ONSPHERE to specify that the quadilateral are on a sphere and 
     23; @keyword DOUBLE use double precision to perform the computation 
     24; 
     25; @keyword ONSPHERE to specify that the quadilateral are on a sphere and 
    2526;    that teir coordinates are longitude-latitude coordinates. In this 
    2627;    case, est-west periodicity, poles singularity and other pbs 
    2728;    related to longitude-latitude coordinates are managed 
    28 ;    automatically.  
     29;    automatically. 
    2930; 
    3031; @keyword ZOOMRADIUS {default=4} 
    … …  
    3233;    zoomradius degree where we look for the the quadrilateral which 
    3334;    contains the (x,y) point) used for the satellite projection 
    34 ;    when /onsphere is activated.  
    35 ;    4 seems to be the minimum which can be used.  
     35;    when /ONSPHERE is activated. 
     36;    4 seems to be the minimum which can be used. 
    3637;    Can be increase if the cell size is larger than 5 degrees. 
    37 ;    
    38 ; @keyword /NOPRINT to suppress the print messages. 
     38; 
     39; @keyword NOPRINT to suppress the print messages. 
    3940; 
    4041; @keyword NEWCOORD 
    … …  
    4546;    quadrilateral number j with (0 <= j <= n_elements(x0)-1) 
    4647; 
    47 ; @restrictions I think degenerated quadrilateral (e.g. flat of 
    48 ; twisted) is not work. This has to be tested. 
    49 ; 
    50 ; @examples  
     48; @restrictions 
     49; I think degenerated quadrilateral (e.g. flat of twisted) is not work. 
     50; This has to be tested. 
     51; 
     52; @examples 
    5153; 
    5254; IDL> x = 1.*[1, 2, 6, 7, 3] 
    … …  
    7072;      Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    7173;      August 2003 
    72 ;      Based on Convert_clic_ij.pro written by Gurvan Madec  
     74;      Based on Convert_clic_ij.pro written by Gurvan Madec 
    7375; 
    7476; @version $Id$ 
    … …  
    100102    x3 = x3 MOD 360 
    101103    x4 = x4 MOD 360 
    102 ; save !map  
     104; save !map 
    103105    save = {map:!map, x:!x, y:!y, z:!z, p:!p} 
    104 ; do a satellite projection  
     106; do a satellite projection 
    105107    IF NOT keyword_set(zoomradius) THEN zoomradius = 4 
    106108    map_set, y[0], x[0], 0, /satellite, sat_p = [1+zoomradius*20/6371.229, 0, 0], /noerase, /iso, /noborder 
    107109; use normal coordinates to reject cells which are out of the projection. 
    108     tmp  = convert_coord(x, y, /DATA, /TO_NORMAL, DOUBLE = double)  
    109     tmp1 = convert_coord(x1, y1, /DATA, /TO_NORMAL, DOUBLE = double)  
    110     tmp2 = convert_coord(x2, y2, /DATA, /TO_NORMAL, DOUBLE = double)  
    111     tmp3 = convert_coord(x3, y3, /DATA, /TO_NORMAL, DOUBLE = double)  
    112     tmp4 = convert_coord(x4, y4, /DATA, /TO_NORMAL, DOUBLE = double)  
     110    tmp  = convert_coord(x, y, /DATA, /TO_NORMAL, DOUBLE = double) 
     111    tmp1 = convert_coord(x1, y1, /DATA, /TO_NORMAL, DOUBLE = double) 
     112    tmp2 = convert_coord(x2, y2, /DATA, /TO_NORMAL, DOUBLE = double) 
     113    tmp3 = convert_coord(x3, y3, /DATA, /TO_NORMAL, DOUBLE = double) 
     114    tmp4 = convert_coord(x4, y4, /DATA, /TO_NORMAL, DOUBLE = double) 
    113115; remove cell which have one corner with coordinates equal to NaN 
    114116    test = finite(tmp1[0, *]+tmp1[1, *]+tmp2[0, *]+tmp2[1, *] $ 
    … …  
    150152; 
    151153    tmp1 = -1 & tmp2 = -1 & tmp3 = -1 & tmp4 = -1 
    152 ; remove cells which are obviously bad  
     154; remove cells which are obviously bad 
    153155    test = (x1 GT x AND x2 GT x AND x3 GT x AND x4 GT x) $ 
    154156      OR (x1 LT x AND x2 LT x AND x3 LT x AND x4 LT x) $ 
    … …  
    179181    ENDIF 
    180182; 
    181     nquad = n_elements(good2)  
     183    nquad = n_elements(good2) 
    182184    x1 = x1[good2] 
    183185    y1 = y1[good2] 
    … …  
    196198;       *((x-x2)*(y3-y2) GT (x3-x2)*(y-y2)) $ 
    197199;       *((x-x3)*(y4-y3) GT (x4-x3)*(y-y3)) $ 
    198 ;       *((x-x4)*(y1-y4) GE (x1-x4)*(y-y4))  
     200;       *((x-x4)*(y1-y4) GE (x1-x4)*(y-y4)) 
    199201; 
    200202; computation of test without any do loop for ntofind points (x,y) and 
    … …  
    202204; test dimensions are (ntofind, nquad) 
    203205; column i of test corresponds to the intersection of point i with all 
    204 ; quadirlateral.  
    205 ; row j of test corresponds to all the points localized in cell j  
     206; quadirlateral. 
     207; row j of test corresponds to all the points localized in cell j 
    206208  test = $ 
    207209; (x-x1) 
    … …  
    277279        ENDIF 
    278280        return,  -1 
    279       END  
     281      END 
    280282      n_elements(found) GT ntofind:BEGIN 
    281283        IF NOT keyword_set(noprint) THEN print, 'The point is in more than one cell' 
    282       END  
     284      END 
    283285      ELSE: 
    284286    ENDCASE 
    285287  ENDIF ELSE BEGIN 
    286 ; if ntofind GT 1, found must be sorted  
     288; if ntofind GT 1, found must be sorted 
    287289; i position of found. this corresponds to one x,y point 
    288290    forsort = found MOD ntofind 

  • trunk/SRC/Interpolation/inrecgrid.pro

    r121 r125  
    11;+ 
    22; 
    3 ; @file_comments given - a list of points, (x,y) position   
    4 ;                - the x and y limits of a rectangular grid 
    5 ;          find in which cell is located each given point. 
     3; @file_comments  
     4, given - a list of points, (x,y) position 
     5;       - the x and y limits of a rectangular grid 
     6; find in which cell is located each given point. 
    67; 
    78; @categories no DO loop, use the wonderfull value_locate function! 
    89; 
    9 ; @param x1d {in}{required}  a 1d array, the x position on the points 
    10 ; @param y1d {in}{required}  a 1d array, the y position on the points 
    11 ; @param left {in}{required} a 1d, monotonically increasing array,  
     10; @param x1d {in}{required} 
     11; a 1d array, the x position on the points 
     12; 
     13; @param y1d {in}{required} 
     14; a 1d array, the y position on the points 
     15; 
     16; @param left {in}{required} 
     17; a 1d, monotonically increasing array, 
    1218; the position of the "left" border of each cell. 
    13 ; @param bottom {in}{required}  a 1d, monotonically increasing array,  
     19; 
     20; @param bottom {in}{required} 
     21; a 1d, monotonically increasing array, 
    1422; the position of the "bottom" border of each cell. 
    1523; 
    16 ; @keyword /output2d to get the output as a 2d array (2,n_elements(x1d)), 
     24; @keyword OUTPUT2D 
     25; to get the output as a 2d array (2,n_elements(x1d)), 
    1726;    with res[0,*] the x index accoring to the 1d array defined by 
    1827;    left and res[1,*] the y index accoring to the 1d array defined by 
    1928;    bottom. 
    2029; 
    21 ; @keyword checkout=[rbgrid,ubgrid] specify the right and upper bondaries of 
     30; @keyword CHECKOUT 
     31; = [rbgrid,ubgrid] specify the right and upper boundaries of 
    2232;    the grid and check if some points are out. 
    2333; 
    24 ; @returns the index on the cell accoring to the 2d array defined by 
    25 ; left and bottom. 
     34; @returns 
     35; the index on the cell accoring to the 2d array defined by left and bottom. 
    2636; 
    27 ; @examples  
     37; @examples 
    2838; 
    29 ;  IDL> a=indgen(5) 
    30 ;  IDL> b=indgen(7) 
    31 ;  IDL> r=inrecgrid([0.25,3.25,2],[4.25,2.8,1.4],a,b) 
    32 ;  IDL> print, r 
     39; IDL> a=indgen(5) 
     40; IDL> b=indgen(7) 
     41; IDL> r=inrecgrid([0.25,3.25,2],[4.25,2.8,1.4],a,b) 
     42; IDL> print, r 
    3343;            20          13           7 
    34 ;  IDL> r=inrecgrid([0.25,3.25,2],[4.25,2.8,1.4],a,a+1,b,b+1,/output2d) 
    35 ;  IDL> print, r 
     44; IDL> r=inrecgrid([0.25,3.25,2],[4.25,2.8,1.4],a,a+1,b,b+1,/output2d) 
     45; IDL> print, r 
    3646;        0.00000      4.00000 
    3747;        3.00000      2.00000 
    3848;        2.00000      1.00000 
    39 ;   
     49; 
    4050; @history 
    41 ;            S. Masson (smasson\@lodyc.jussieu.fr) 
    42 ;                      July 3rd, 2002 
    43 ;                      October 3rd, 2003: use value_locate 
     51; S. Masson (smasson\@lodyc.jussieu.fr) 
     52; July 3rd, 2002 
     53; October 3rd, 2003: use value_locate 
    4454; 
    4555; @version $Id$ 
    4656; 
    4757;- 
    48  
    49 FUNCTION inrecgrid, x1d, y1d, left, bottom, output2d = output2d, checkout = checkout 
     58FUNCTION inrecgrid, x1d, y1d, left, bottom, OUTPUT2D = output2d, CHECKOUT = checkout 
    5059; 
    5160  compile_opt idl2, strictarrsubs 
    … …  
    7180  out = where(xpos EQ -1 OR ypos EQ -1) 
    7281  IF out[0] NE -1 THEN res[out] = -1 
    73 ;   
     82; 
    7483  RETURN, res 
    7584 

  • trunk/SRC/Interpolation/ll_narcs_distances.pro

    r121 r125  
    33; @file_comments 
    44; This function returns the longitude and latitude [lon, lat] of 
    5 ;a point a given arc distance (-pi <= Arc_Dist <= pi), and azimuth (Az), 
    6 ;from a specified location Lon0, lat0. 
    7 ;       Same as LL_ARC_DISTANCE but for n points without do loop. 
     5; a point a given arc distance (-pi <= Arc_Dist <= pi), and azimuth (Az), 
     6; from a specified location Lon0, lat0. 
     7; Same as LL_ARC_DISTANCE but for n points without do loop. 
     8; 
     9; Formula from Map Projections - a working manual.  USGS paper 
     10; 1395.  Equations (5-5) and (5-6). 
    811; 
    912; @categories Mapping, geography 
    1013; 
    11 ;    @param Lon0 {in}{required}  An array containing the longitude of the starting point. 
    12 ;             Values are assumed to be in radians unless the keyword 
    13 ;             DEGREES is set. 
    14 ;    @param Lat0 {in}{required}  An array containing the latitude of the starting point. 
    15 ;             Values are assumed to be in radians unless the keyword 
    16 ;             DEGREES is set. 
    17 ;    @param Arc_Dist {in}{required}  The arc distance from Lon_lat0. The value must be between 
     14; @param Lon0 {in}{required} 
     15; An array containing the longitude of the starting point. 
     16; Values are assumed to be in radians unless the keyword DEGREES is set. 
     17; 
     18; @param Lat0 {in}{required} 
     19; An array containing the latitude of the starting point. 
     20; Values are assumed to be in radians unless the keyword DEGREES is set. 
     21; 
     22; @param Arc_Dist {in}{required} 
     23; The arc distance from Lon_lat0. The value must be between 
    1824; -!PI and +!PI. To express distances in arc units, divide 
    1925;  by the radius of the globe expressed in the original units. 
    2026;  For example, if the radius of the earth is 6371 km, divide 
    21 ;  the distance in km by 6371 to obtain the arc distance.     
    22 ;    @param Az {in}{required}   The azimuth from Lon_lat0. The value is assumed to be in 
    23 ;  radians unless the keyword DEGREES is set. 
     27;  the distance in km by 6371 to obtain the arc distance. 
    2428; 
    25 ; @keyword    DEGREES  Set this keyword to express all measurements and 
    26 ;  results in degrees. 
     29; @param Az {in}{required} 
     30; The azimuth from Lon_lat0. The value is assumed to be in 
     31; radians unless the keyword DEGREES is set. 
     32; 
     33; @keyword DEGREES 
     34; Set this keyword to express all measurements and results in degrees. 
    2735; 
    2836; @returns 
    29 ; a (2, n) array containing the  
    30 ;       longitude / latitude of the resultings points. Values are in radians 
    31 ;       unless the keyword DEGREES is set. 
     37; a (2, n) array containing the longitude/latitude of the resultings points. 
     38; Values are in radians unless the keyword DEGREES is set. 
    3239; 
    33 ; @file_comments 
    34 ;Formula from Map Projections - a working manual.  USGS paper 
    35 ;1395.  Equations (5-5) and (5-6). 
    36 ; 
    37 ; @examples  
     40; @examples 
    3841; IDL> Lon_lat0 = [1.0, 2.0]; Initial point specified in radians 
    3942; IDL> Arc_Dist = 2.0; Arc distance in radians 
    … …  
    4346;       2.91415    -0.622234 
    4447; 
    45 ;IDL> lon0 = [-10, 20, 100] 
    46 ;IDL> lat0 = [0, -10, 45] 
    47 ;IDL> lon1 = [10, 60, 280] 
    48 ;IDL> lat1 = [0, 10, 45] 
    49 ;IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi, /two_by_two) 
    50 ;IDL> earthradius = 6378206.4d0 
    51 ;IDL> res = ll_narcs_distances(lon0, lat0, dist/earthradius, azi, /degrees) 
    52 ;IDL> print, reform(res[0, *]) 
     48; IDL> lon0 = [-10, 20, 100] 
     49; IDL> lat0 = [0, -10, 45] 
     50; IDL> lon1 = [10, 60, 280] 
     51; IDL> lat1 = [0, 10, 45] 
     52; IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi, /two_by_two) 
     53; IDL> earthradius = 6378206.4d0 
     54; IDL> res = ll_narcs_distances(lon0, lat0, dist/earthradius, azi, /degrees) 
     55; IDL> print, reform(res[0, *]) 
    5356;       10.000000       60.000000       280.00000 
    54 ;IDL> print, reform(res[1, *]) 
    55 ;          1.1999280e-15       10.000000       45.000000 
     57; IDL> print, reform(res[1, *]) 
     58;           1.1999280e-15       10.000000       45.000000 
    5659; 
    5760; @history 
    … …  
    6265; @version $Id$ 
    6366; 
    64 ;- 
    65  
    66 ;+ 
    67 ; @file_comments Return the [lon, lat] of the point a given arc distance  
    68 ;(-pi <= arc_dist <= pi), 
    69 ; and azimuth (az), from lon_lat0. 
    7067;- 
    7168; 
    … …  
    9491 
    9592  zero = where(arc_dist eq 0, count) 
    96   IF count NE 0 THEN BEGIN  
     93  IF count NE 0 THEN BEGIN 
    9794    lam[zero] = lon0[zero] 
    9895    phi[zero] = lat0[zero] 
    99   ENDIF  
     96  ENDIF 
    10097 
    10198  if keyword_set(degs) then return, transpose([[lam], [phi]]) / s $ 
    … …  
    103100 
    104101end 
    105  
    106  

  • trunk/SRC/Interpolation/map_npoints.pro

    r121 r125  
    1 ;+ 
    2 ; 
    3 ; @file_comments 
    4 ;Return the distance in meter between all np0 points P0 and all 
    5 ;       np1 points P1 on a sphere. If keyword /TWO_BY_TWO is given then 
    6 ;       returns the distances between number n of P0 points and number 
    7 ;       n of P1 points (in that case, np0 and np1 must be equal). 
    8 ;       Same as map_2points with the meter parameter but for n points 
    9 ;       without do loop. 
     1;+  
     2;  
     3; @file_comments  
     4; Return the distance in meter between all np0 points P0 and all  
     5; np1 points P1 on a sphere. If keyword /TWO_BY_TWO is given then  
     6; returns the distances between number n of P0 points and number  
     7; n of P1 points (in that case, np0 and np1 must be equal). 
     8; Same as map_2points with the meter parameter but for n points 
     9; without do loop. 
    1010; 
    1111; @categories Maps 
    1212; 
    1313; @param Lon0 {in}{required} 
    14 ; @param Lat0 {in}{required}  
    15 ; np0 elements vector. longitudes and latitudes of np0 points P0  
     14; @param Lat0 {in}{required} 
     15; np0 elements vector. longitudes and latitudes of np0 points P0 
    1616; 
    1717; @param Lon1 {in}{required} 
    18 ; @param Lat1 {in}{required}   
    19 ; np1 elements vector. longitude and latitude of np1 points P1  
     18; @param Lat1 {in}{required} 
     19; np1 elements vector. longitude and latitude of np1 points P1 
    2020; 
    21 ; @keyword AZIMUTH A named variable that will receive the azimuth of the great 
    22 ;       circle  connecting the two points, P0 to P1 
    23 ; @keyword /MIDDLE to get the longitude/latitude of the middle point betwen P0 and P1. 
    24 ; @keyword RADIANS if set, inputs and angular outputs are in radians, otherwise 
    25 ;degrees. 
    26 ; @keyword RADIUS {default=6378206.4d0}  
    27 ; If given, return the distance between the two points calculated using the  
     21; @keyword AZIMUTH 
     22; A named variable that will receive the azimuth of the great 
     23; circle connecting the two points, P0 to P1 
     24; 
     25; @keyword MIDDLE 
     26; to get the longitude/latitude of the middle point betwen P0 and P1. 
     27; 
     28; @keyword RADIANS 
     29; if set, inputs and angular outputs are in radians, otherwise degrees. 
     30; 
     31; @keyword RADIUS {default=6378206.4d0} 
     32; If given, return the distance between the two points calculated using the 
    2833; given radius. 
    2934; Default value is the Earth radius. 
    3035; 
    31 ; @keyword TWO_BY_TWO If given,then Map_nPoints returns the distances between 
    32 ;       number n of P0 points and number n of P1 points (in that case, 
    33 ;       np0 and np1 must be equal). 
     36; @keyword TWO_BY_TWO 
     37; If given,then Map_nPoints returns the distances between number n of 
     38; P0 points and number n of P1 points 
     39; In that case, np0 and np1 must be equal. 
    3440; 
    3541; @returns 
    36 ;       An (np0,np1) array giving the distance in meter between np0 
    37 ;       points P0 and np1 points P1. Element (i,j) of the ouput is the 
    38 ;       distance between element P0[i] and P1[j]. 
    39 ;       If keyword /TWO_BY_TWO is given then Map_nPoints returns 
    40 ;       an np-element vector giving the distance in meter between P0[i] 
    41 ;       and P1[i] (in that case, we have np0 = np1 = np) 
    42 ;       if /MIDDLE see this keyword. 
     42; An (np0,np1) array giving the distance in meter between np0 
     43; points P0 and np1 points P1. Element (i,j) of the ouput is the 
     44; distance between element P0[i] and P1[j]. 
     45; If keyword /TWO_BY_TWO is given then Map_nPoints returns 
     46; an np-element vector giving the distance in meter between P0[i] 
     47; and P1[i] (in that case, we have np0 = np1 = np) ; if /MIDDLE see this keyword. 
    4348; 
    4449; @examples 
    4550; IDL> print, $ 
    46 ; map_npoints([-105.15,1],[40.02,1],[-0.07,100,50],[51.30,20,0]) 
    47 ;       7551369.3      5600334.8 
    48 ;       12864354.      10921254. 
    49 ;       14919237.      5455558.8 
     51; IDL> map_npoints([-105.15,1],[40.02,1],[-0.07,100,50],[51.30,20,0]) 
     52; 7551369.3 5600334.8 
     53; 12864354. 10921254. 
     54; 14919237. 5455558.8 
    5055; 
    5156; IDL> lon0 = [-10, 20, 100] 
    … …  
    5358; IDL> lon1 = [10, 60, 280] 
    5459; IDL> lat1 = [0, 10, 45] 
    55 ; IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi) 
     60; IDL> dist = map_npoints(lon0, lat0, lon1, lat1, AZIMUTH = azi) 
    5661; IDL> help, dist, azi 
    57 ; DIST            DOUBLE    = Array[3, 3] 
    58 ; AZI             DOUBLE    = Array[3, 3] 
     62; DIST DOUBLE = Array[3, 3] 
     63; AZI DOUBLE = Array[3, 3] 
    5964; IDL> print, dist[4*lindgen(3)], azi[4*lindgen(3)] 
    60 ;       2226414.0       4957944.5      10018863. 
    61 ;       90.000000       64.494450  4.9615627e-15 
    62 ; IDL> dist = map_npoints(lon0, lat0, lon1, lat1, azimuth = azi, /two_by_two) 
     65; 2226414.0 4957944.5 10018863. 
     66; 90.000000 64.494450 4.9615627e-15 
     67; IDL> dist = map_npoints(lon0, lat0, lon1, lat1, AZIMUTH = azi, /TWO_BY_TWO) 
    6368; IDL> help, dist, azi 
    64 ; DIST            DOUBLE    = Array[3] 
    65 ; AZI             DOUBLE    = Array[3] 
     69; DIST DOUBLE = Array[3] 
     70; AZI DOUBLE = Array[3] 
    6671; IDL> print, dist, azi 
    67 ;       2226414.0       4957944.5      10018863. 
    68 ;       90.000000       64.494450  4.9615627e-15 
     72; 2226414.0 4957944.5 10018863. 
     73; 90.000000 64.494450 4.9615627e-15 
    6974; IDL> print, map_2points(lon0[0], lat0[0], lon1[0], lat1[0]) 
    70 ;       20.000000      90.000000 
    71 ; IDL> print, map_npoints(lon0[0], lat0[0], lon1[0], lat1[0], azi=azi)/6378206.4d0 / !dtor, azi 
    72 ;       20.000000 
    73 ;       90.000000 
     75; 20.000000 90.000000 
     76; IDL> print, map_npoints(lon0[0], lat0[0], lon1[0], lat1[0], AZIMUTH=azi)/6378206.4d0 / !dtor, azi 
     77; 20.000000 
     78; 90.000000 
    7479; 
    7580; IDL> lon0 = [-10, 20, 100] 
    … …  
    7782; IDL> lon1 = [10, 60, 280] 
    7883; IDL> lat1 = [0, 10, 45] 
    79 ; IDL> mid = map_npoints(lon0, lat0, lon1, lat1, /middle, /two_by_two) 
     84; IDL> mid = map_npoints(lon0, lat0, lon1, lat1, /MIDDLE, /TWO_BY_TWO) 
    8085; IDL> print, reform(mid[0,*]), reform(mid[1,*]) 
    81 ;       0.0000000       40.000000      190.00000 
    82 ;       0.0000000  -1.5902773e-15      90.000000 
     86; 0.0000000 40.000000 190.00000 
     87; 0.0000000 -1.5902773e-15 90.000000 
    8388; IDL> print, (map_2points(lon0[0], lat0[0], lon1[0], lat1[0], npath = 3))[*, 1] 
    84 ;       0.0000000      0.0000000 
     89; 0.0000000 0.0000000 
    8590; IDL> print, (map_2points(lon0[1], lat0[1], lon1[1], lat1[1], npath = 3))[*, 1] 
    86 ;       40.000000 -1.5902773e-15 
     91; 40.000000 -1.5902773e-15 
    8792; IDL> print, (map_2points(lon0[2], lat0[2], lon1[2], lat1[2], npath = 3))[*, 1] 
    88 ;       190.00000      90.000000 
     93; 190.00000 90.000000 
    8994; 
    9095; @history 
    91 ;       Based on the IDL function map_2points.pro,v 1.6 2001/01/15 
     96; Based on the IDL function map_2points.pro,v 1.6 2001/01/15 
    9297; Sebastien Masson (smasson\@lodyc.jussieu.fr) 
    93 ;                  October 2003 
     98; October 2003 
    9499; 
    95100; @version $Id$ 
    96101; 
    97102;- 
    98 Function Map_npoints, lon0, lat0, lon1, lat1, azimuth = azimuth $ 
    99   , RADIANS = radians, RADIUS = radius, MIDDLE = middle, TWO_BY_TWO = two_by_two 
     103Function Map_npoints, lon0, lat0, lon1, lat1, AZIMUTH = azimuth $ 
     104 , RADIANS = radians, RADIUS = radius, MIDDLE = middle, TWO_BY_TWO = two_by_two 
    100105 
    101   COMPILE_OPT idl2, strictarrsubs 
     106 COMPILE_OPT idl2, strictarrsubs 
    102107 
    103   IF (N_PARAMS() LT 4) THEN $ 
    104     MESSAGE, 'Incorrect number of arguments.' 
     108 IF (N_PARAMS() LT 4) THEN $ 
     109 MESSAGE, 'Incorrect number of arguments.' 
    105110 
    106   np0 = n_elements(lon0)  
    107   IF n_elements(lat0) NE np0 THEN $ 
    108     MESSAGE, 'lon0 and lat0 must have the same number of elements' 
    109   np1 = n_elements(lon1)  
    110   IF n_elements(lat1) NE np1 THEN $ 
    111     MESSAGE, 'lon1 and lat1 must have the same number of elements' 
    112   if keyword_set(two_by_two) AND np0 NE np1 then $ 
    113     MESSAGE, 'When using two_by_two keyword, P0 and P1 must have the same number of elements' 
     111 np0 = n_elements(lon0) 
     112 IF n_elements(lat0) NE np0 THEN $ 
     113 MESSAGE, 'lon0 and lat0 must have the same number of elements' 
     114 np1 = n_elements(lon1) 
     115 IF n_elements(lat1) NE np1 THEN $ 
     116 MESSAGE, 'lon1 and lat1 must have the same number of elements' 
     117 if keyword_set(two_by_two) AND np0 NE np1 then $ 
     118 MESSAGE, 'When using two_by_two keyword, P0 and P1 must have the same number of elements' 
    114119 
    115   mx = MAX(ABS([lat0[*], lat1[*]])) 
    116   pi2 = !dpi/2 
    117   IF (mx GT (KEYWORD_SET(radians) ? pi2 : 90)) THEN $ 
    118     MESSAGE, 'Value of Latitude is out of allowed range.' 
     120 mx = MAX(ABS([lat0[*], lat1[*]])) 
     121 pi2 = !dpi/2 
     122 IF (mx GT (KEYWORD_SET(radians) ? pi2 : 90)) THEN $ 
     123 MESSAGE, 'Value of Latitude is out of allowed range.' 
    119124 
    120   k = KEYWORD_SET(radians) ? 1.0d0 : !dpi/180.0 
     125 k = KEYWORD_SET(radians) ? 1.0d0 : !dpi/180.0 
    121126;Earth equatorial radius, meters, Clarke 1866 ellipsoid 
    122   r_sphere =  n_elements(RADIUS) NE 0 ? RADIUS : 6378206.4d0  
     127 r_sphere = n_elements(RADIUS) NE 0 ? RADIUS : 6378206.4d0 
    123128; 
    124   coslt1 = cos(k*lat1[*]) 
    125   sinlt1 = sin(k*lat1[*]) 
    126   coslt0 = cos(k*lat0[*]) 
    127   sinlt0 = sin(k*lat0[*]) 
     129 coslt1 = cos(k*lat1[*]) 
     130 sinlt1 = sin(k*lat1[*]) 
     131 coslt0 = cos(k*lat0[*]) 
     132 sinlt0 = sin(k*lat0[*]) 
    128133; 
    129   IF np0 EQ np1 AND np1 EQ 1 THEN two_by_two = 1 
     134 IF np0 EQ np1 AND np1 EQ 1 THEN two_by_two = 1 
    130135; 
    131   if NOT keyword_set(two_by_two) THEN BEGIN  
    132     coslt1 = replicate(1.0d0, np0)#temporary(coslt1) 
    133     sinlt1 = replicate(1.0d0, np0)#temporary(sinlt1) 
    134     coslt0 = temporary(coslt0)#replicate(1.0d0, np1) 
    135     sinlt0 = temporary(sinlt0)#replicate(1.0d0, np1) 
    136   ENDIF  
     136 if NOT keyword_set(two_by_two) THEN BEGIN 
     137 coslt1 = replicate(1.0d0, np0)#temporary(coslt1) 
     138 sinlt1 = replicate(1.0d0, np0)#temporary(sinlt1) 
     139 coslt0 = temporary(coslt0)#replicate(1.0d0, np1) 
     140 sinlt0 = temporary(sinlt0)#replicate(1.0d0, np1) 
     141 ENDIF 
    137142; 
    138   if keyword_set(two_by_two) THEN BEGIN  
    139     cosl0l1 = cos(k*(lon1[*]-lon0[*])) 
    140     sinl0l1 = sin(k*(lon1[*]-lon0[*])) 
    141   ENDIF ELSE BEGIN  
    142     cosl0l1 = cos(k*(replicate(1.0d0, np0)#lon1[*]-lon0[*]#replicate(1.0d0, np1))) 
    143     sinl0l1 = sin(k*(replicate(1.0d0, np0)#lon1[*]-lon0[*]#replicate(1.0d0, np1))) 
    144   ENDELSE  
     143 if keyword_set(two_by_two) THEN BEGIN 
     144 cosl0l1 = cos(k*(lon1[*]-lon0[*])) 
     145 sinl0l1 = sin(k*(lon1[*]-lon0[*])) 
     146 ENDIF ELSE BEGIN 
     147 cosl0l1 = cos(k*(replicate(1.0d0, np0)#lon1[*]-lon0[*]#replicate(1.0d0, np1))) 
     148 sinl0l1 = sin(k*(replicate(1.0d0, np0)#lon1[*]-lon0[*]#replicate(1.0d0, np1))) 
     149 ENDELSE 
    145150 
    146   cosc = sinlt0 * sinlt1 + coslt0 * coslt1 * cosl0l1 ;Cos of angle between pnts 
     151 cosc = sinlt0 * sinlt1 + coslt0 * coslt1 * cosl0l1 ;Cos of angle between pnts 
    147152; Avoid roundoff problems by clamping cosine range to [-1,1]. 
    148   cosc = -1.0d0 > cosc < 1.0d0 
     153 cosc = -1.0d0 > cosc < 1.0d0 
    149154; 
    150   if arg_present(azimuth) OR keyword_set(middle) then begin 
    151     sinc = sqrt(1.0d0 - cosc*cosc) 
    152     bad = where(abs(sinc) le 1.0e-7) 
    153     IF bad[0] NE -1 THEN sinc[bad] = 1 
    154     cosaz = (coslt0 * sinlt1 - sinlt0*coslt1*cosl0l1) / sinc  
    155     sinaz = sinl0l1*coslt1/sinc 
    156     IF bad[0] NE -1 THEN BEGIN  
    157       sinc[bad] = 0.0d0 
    158       sinaz[bad] = 0.0d0 
    159       cosaz[bad] = 1.0d0 
    160     ENDIF 
    161   ENDIF 
     155 if arg_present(azimuth) OR keyword_set(middle) then begin 
     156 sinc = sqrt(1.0d0 - cosc*cosc) 
     157 bad = where(abs(sinc) le 1.0e-7) 
     158 IF bad[0] NE -1 THEN sinc[bad] = 1 
     159 cosaz = (coslt0 * sinlt1 - sinlt0*coslt1*cosl0l1) / sinc 
     160 sinaz = sinl0l1*coslt1/sinc 
     161 IF bad[0] NE -1 THEN BEGIN 
     162 sinc[bad] = 0.0d0 
     163 sinaz[bad] = 0.0d0 
     164 cosaz[bad] = 1.0d0 
     165 ENDIF 
     166 ENDIF 
    162167; 
    163   IF keyword_set(middle) then BEGIN 
     168 IF keyword_set(middle) then BEGIN 
    164169 
    165     s0 = 0.5d0 * acos(cosc) 
     170 s0 = 0.5d0 * acos(cosc) 
    166171 ; 
    167     coss = cos(s0) 
    168     sins = sin(s0) 
    169 ;     
    170     lats = asin(sinlt0 * coss + coslt0 * sins * cosaz) / k 
    171     lons = atan(sins * sinaz, coslt0 * coss - sinlt0 * sins * cosaz) / k 
     172 coss = cos(s0) 
     173 sins = sin(s0) 
    172174; 
    173     if keyword_set(two_by_two) THEN BEGIN  
    174       return, transpose([[lon0[*] + lons], [lats]]) 
    175     ENDIF ELSE BEGIN  
    176       return, [ [[lon0[*]#replicate(1.0d0, np1) + lons]], [[lats]] ] 
    177     ENDELSE  
     175 lats = asin(sinlt0 * coss + coslt0 * sins * cosaz) / k 
     176 lons = atan(sins * sinaz, coslt0 * coss - sinlt0 * sins * cosaz) / k 
    178177; 
    179   ENDIF 
     178 if keyword_set(two_by_two) THEN BEGIN 
     179 return, transpose([[lon0[*] + lons], [lats]]) 
     180 ENDIF ELSE BEGIN 
     181 return, [ [[lon0[*]#replicate(1.0d0, np1) + lons]], [[lats]] ] 
     182 ENDELSE 
    180183; 
    181   if arg_present(azimuth) then begin 
    182     azimuth = atan(sinaz, cosaz) 
    183     IF k NE 1.0d0 THEN azimuth = temporary(azimuth) / k 
    184    ENDIF 
     184 ENDIF 
     185; 
     186 if arg_present(azimuth) then begin 
     187 azimuth = atan(sinaz, cosaz) 
     188 IF k NE 1.0d0 THEN azimuth = temporary(azimuth) / k 
     189 ENDIF 
    185190 return, acos(cosc) * r_sphere 
    186191; 

  • trunk/SRC/Interpolation/neighbor.pro

    r121 r125  
    22; 
    33; @file_comments 
    4 ;find the closetest point of (P0) within a list of np1 points 
    5 ;P1 Which can be on a sphere  
     4; find the closetest point of (P0) within a list of np1 points 
     5; P1 Which can be on a sphere 
    66; 
    77; @categories Maps 
    88; 
    9 ; @param p0lon {in}{required}  scalar. longitudes of point P0.  
    10 ; @param p0lat {in}{required}  scalar. latitudes of point P0.  
    11 ; @param neighlon {in}{optional}  
    12 ; @param neighlat {in}{optional}  
     9; @param p0lon {in}{required}  scalar. longitudes of point P0. 
     10; @param p0lat {in}{required}  scalar. latitudes of point P0. 
     11; @param neighlon {in}{optional} 
     12; @param neighlat {in}{optional} 
    1313; 
    14 ; @keyword RADIANS if set, inputs and angular outputs are in radians, otherwise 
    15 ;degrees. 
    16 ; @keyword DISTANCE dis, to get back the distances between P0 and the np1 
    17 ;   points P1 in the variable dis. 
    18 ; @keyword /SPHERE to activate if points are located on a sphere. 
     14; @keyword RADIANS 
     15; if set, inputs and angular outputs are in radians, otherwise degrees. 
     16; 
     17; @keyword DISTANCE 
     18; dis, to get back the distances between P0 and the np1 points P1 in the 
     19; variable dis. 
     20; 
     21; @keyword SPHERE to activate if points are located on a sphere. 
    1922; 
    2023; @returns 
    21 ;       index giving the P1[index] point that is the closest point of (P0) 
     24; index giving the P1[index] point that is the closest point of (P0) 
    2225; 
    2326; @examples 
    24 ;       IDL> print, neighbor(-105.15,40.02,[-0.07,100,50],[51.30,20,0], $ 
    25 ;           distance=dis) 
     27; IDL> print, neighbor(-105.15,40.02,[-0.07,100,50],[51.30,20,0], $ 
     28; IDL> distance=dis) 
    2629;                  0 
    27 ;       IDL> print, dis 
     30; IDL> print, dis 
    2831;             105.684      206.125      160.228 
    2932; 
    … …  
    4447  IF  n_elements(p0lat) NE 1 THEN MESSAGE, 'Sorry p0lat must be a scalar' 
    4548  p0lat = p0lat[0] 
    46   nneig = n_elements(neighlon)  
     49  nneig = n_elements(neighlon) 
    4750  IF  n_elements(neighlat) NE nneig  THEN $ 
    4851    MESSAGE, 'neighlon and neighlat must have the same number of elements' 
    … …  
    5255    distance = Map_nPoints(p0lon, p0lat, neighlon, neighlat $ 
    5356                       , radius = radius, radians = radians) 
    54   ENDIF ELSE BEGIN  
     57  ENDIF ELSE BEGIN 
    5558    distance = (neighlon-p0lon)^2+(neighlat-p0lat)^2 
    5659    IF arg_present(distance) THEN distance = sqrt(distance) 

  • trunk/SRC/Interpolation/quadrilateral2square.pro

    r121 r125  
    11;+ 
    22; 
    3 ; @file_comments warm (or map) an arbitrary quadrilateral onto a unit square  
     3; @file_comments 
     4; warm (or map) an arbitrary quadrilateral onto a unit square 
    45; according to the 4-point correspondences: 
    56;       (x0,y0) -> (0,0) 
    … …  
    3536; 
    3637;     (2,n) array: the new coodinates (xout, yout) of the (xin,yin) 
    37 ;     point(s) after mapping.  
     38;     point(s) after mapping. 
    3839;     If xin is a scalar, then n is equal to the number of elements of 
    3940;     x0. If xin is an array , then n is equal to the number of 
    4041;     elements of xin. 
    4142; 
    42 ; @restrictions I think degenerated quadrilateral (e.g. flat of 
    43 ; twisted) is not work. This has to be tested. 
     43; @restrictions 
     44; I think degenerated quadrilateral (e.g. flat of twisted) is not work. 
     45; This has to be tested. 
    4446; 
    45 ; @examples  
     47; @examples 
    4648; 
    4749; IDL> splot,[0,5],[0,3],/nodata,xstyle=1,ystyle=1 
    … …  
    6163;      IEEE Computer Society Press, Los Alamitos, California 
    6264;      Chapter 3, see p 52-56 
    63 ;       
     65; 
    6466; 
    6567; @version $Id$ 
    … …  
    109111; 
    110112; compute the adjoint matrix 
    111 ;  
     113; 
    112114  IF keyword_set(double) THEN adj = dblarr(9, n_elements(x0)) $ 
    113115  ELSE adj = fltarr(9, n_elements(x0)) 
    … …  
    122124  adj[7, *] = a[6, *]*a[1, *]-a[0, *]*a[7, *] 
    123125  adj[8, *] = a[0, *]*a[4, *]-a[3, *]*a[1, *] 
    124 ;   
     126; 
    125127  IF n_elements(xin) EQ 1 THEN BEGIN 
    126     xin = replicate(xin, n_elements(x0))  
    127     yin = replicate(yin, n_elements(x0))  
     128    xin = replicate(xin, n_elements(x0)) 
     129    yin = replicate(yin, n_elements(x0)) 
    128130  ENDIF 
    129131; 

  • trunk/SRC/Interpolation/spl_fstdrv.pro

    r121 r125  
    44;+ 
    55; 
    6 ; @file_comments SPL_FSTDRV returns the values of the first derivative of 
     6; @file_comments 
     7; SPL_FSTDRV returns the values of the first derivative of 
    78; the interpolating function at the points X2i. it is a double 
    89; precision array. 
    … …  
    1415; in a way that interpolated value are also in ascending order 
    1516; 
    16 ; @examples  
     17; @examples 
    1718; IDL> y2 =  spl_fstdrv(x, y, yscd, x2) 
    1819; 
    19 ;    @param x {in}{required}  An n-element (at least 2) input vector that specifies the 
    20 ;    tabulate points in ascending order. 
     20; @param x {in}{required} 
     21; An n-element (at least 2) input vector that specifies the 
     22; tabulate points in ascending order. 
    2123; 
    22 ;    @param y {in}{required}  f(x) = y. An n-element input vector that specifies the values 
    23 ;    of the tabulated function F(Xi) corresponding to Xi. 
     24; @param y {in}{required} 
     25; f(x) = y. An n-element input vector that specifies the values 
     26; of the tabulated function F(Xi) corresponding to Xi. 
    2427; 
    25 ;    @param yscd {in}{required}  The output from SPL_INIT for the specified X and Y. 
     28; @param yscd {in}{required} 
     29; The output from SPL_INIT for the specified X and Y. 
    2630; 
    27 ;    @param x2 {in}{required}  The input values for which the first derivative values are 
    28 ;    desired. X can be scalar or an array of values. 
     31; @param x2 {in}{required} 
     32; The input values for which the first derivative values are desired. 
     33; X can be scalar or an array of values. 
    2934 
    30 ; @returns  
     35; @returns 
    3136; 
    32 ;    y2: f'(x2) = y2.  
     37;    y2: f'(x2) = y2. 
    3338; 
    3439; @history 
    … …  
    5055  ny = n_elements(y) 
    5156; x must have at least 2 elements 
    52   IF nx LT 2 THEN stop    
     57  IF nx LT 2 THEN stop 
    5358; y must have the same number of elements than x 
    5459  IF nx NE ny THEN stop 
    55 ; define loc in a way that  
     60; define loc in a way that 
    5661;  if loc[i] eq -1   :                 x2[i] <  x[0] 
    5762;  if loc[i] eq nx2-1:                 x2[i] >= x[nx-1] 
    5863;  else              :    x[loc[i]] <= x2[i] <  x[loc[i]+1] 
    5964  loc = value_locate(x, x2) 
    60 ; change loc in order to  
     65; change loc in order to 
    6166; use x[0]    and x[1]    even if x2[i] <  x[0]    -> extrapolation 
    6267; use x[nx-2] and x[nx-1] even if x2[i] >= x[nx-1] -> extrapolation 
    63   loc = 0 > temporary(loc) < (nx-2)  
     68  loc = 0 > temporary(loc) < (nx-2) 
    6469; distance between to consecutive x 
    6570  deltax = x[loc+1]-x[loc] 
    … …  
    7479    - 1.0d/6.0d * (3.0d*a*a - 1.0d) * deltax * yscd[loc] $ 
    7580    + 1.0d/6.0d * (3.0d*b*b - 1.0d) * deltax * yscd[loc+1] 
    76 ; beware of the computation precision...  
     81; beware of the computation precision... 
    7782; force near zero values to be exactly 0.0 
    7883  zero = where(abs(yfrst) LT 1.e-10) 

  • trunk/SRC/Interpolation/spl_incr.pro

    r121 r125  
    1212; in a way that interpolated values are also monotonically increasing. 
    1313; 
    14 ; @param x1 {in}{required}   
    15 ; An n-element (at least 2) input vector that specifies the tabulate points in  
     14; @param x1 {in}{required} 
     15; An n-element (at least 2) input vector that specifies the tabulate points in 
    1616; a strict ascending order. 
    1717; 
    18 ; @param y1 {in}{required}   
     18; @param y1 {in}{required} 
    1919; f(x) = y. An n-element input vector that specifies the values 
    2020;    of the tabulated function F(Xi) corresponding to Xi. As f is 
    … …  
    2222;    monotonically increasing. y can have equal consecutive values. 
    2323; 
    24 ; @param x2 {in}{required}   
     24; @param x2 {in}{required} 
    2525; The input values for which the interpolated values are 
    26 ; desired. Its values must be strictly monotonically increasing.  
     26; desired. Its values must be strictly monotonically increasing. 
    2727; 
    2828; @param der2 
    29 ; @param x  
    30 ; 
    31 ; @returns  
     29; @param x 
     30; 
     31; @returns 
    3232; 
    3333;    y2: f(x2) = y2. Double precision array 
    … …  
    3737;   values (amplitude smaller than 1.e-6)... 
    3838; 
    39 ; @examples  
     39; @examples 
    4040; 
    4141; IDL> n = 100L 
    42 ; IDL> x = (dindgen(n))^2  
     42; IDL> x = (dindgen(n))^2 
    4343; IDL> y = abs(randomn(0, n)) 
    4444; IDL> y[n/2:n/2+1] = 0. 
    … …  
    5353; IDL> oplot, x2, y2, color = 100 
    5454; IDL> c = y2[1:n2-1] - y2[0:n2-2] 
    55 ; IDL> print, min(c) LT 0  
     55; IDL> print, min(c) LT 0 
    5656; IDL> print, min(c, max = ma), ma 
    5757; IDL> splot,c,xstyle=1,ystyle=1, yrange=[-.01,.05], ysurx=.25, petit = [1, 2, 2], /noerase 
    … …  
    8787 
    8888;+ 
    89 ; @param x1 {in}{required}   
    90 ; An n-element (at least 2) input vector that specifies the tabulate points in  
     89; @param x1 {in}{required} 
     90; An n-element (at least 2) input vector that specifies the tabulate points in 
    9191; a strict ascending order. 
    9292; 
    93 ; @param y1 {in}{required}   
     93; @param y1 {in}{required} 
    9494; f(x) = y. An n-element input vector that specifies the values 
    9595;    of the tabulated function F(Xi) corresponding to Xi. As f is 
    … …  
    9797;    monotonically increasing. y can have equal consecutive values. 
    9898; 
    99 ; @param x2 {in}{required}   
     99; @param x2 {in}{required} 
    100100; The input values for which the interpolated values are 
    101 ; desired. Its values must be strictly monotonically increasing.  
     101; desired. Its values must be strictly monotonically increasing. 
    102102; 
    103103; @param der2 
    104 ; @param x  
     104; @param x 
    105105; 
    106106;- 
    … …  
    134134; @keyword YPN_1 The first derivative of the interpolating function at the 
    135135;    point Xn-1. If YPN_1 is omitted, the second derivative at the 
    136 ;    boundary is set to zero, resulting in a "natural spline."  
     136;    boundary is set to zero, resulting in a "natural spline." 
    137137;- 
    138138FUNCTION spl_incr, x, y, x2, YP0 = yp0, YPN_1 = ypn_1 
    … …  
    148148  nx2 = n_elements(x2) 
    149149; x must have at least 2 elements 
    150   IF nx LT 2 THEN stop  
     150  IF nx LT 2 THEN stop 
    151151; y must have the same number of elements than x 
    152152  IF nx NE ny THEN stop 
    153153; x be monotonically increasing 
    154   IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop  
     154  IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop 
    155155; x2 be monotonically increasing 
    156156  IF N_ELEMENTS(X2) GE 2 THEN $ 
    157   IF min(x2[1:nx2-1]-x2[0:nx2-2])  LE 0 THEN stop  
     157  IF min(x2[1:nx2-1]-x2[0:nx2-2])  LE 0 THEN stop 
    158158; y be monotonically increasing 
    159   IF min(y[1:ny-1]-y[0:ny-2]) LT 0 THEN stop  
     159  IF min(y[1:ny-1]-y[0:ny-2]) LT 0 THEN stop 
    160160;--------------------------------- 
    161161; first check: check if two consecutive values are equal 
    … …  
    172172    xinx2_1 = value_locate(x2, x[bad+1]) 
    173173; 
    174 ; left side ... if there is x2 values smaller that x[bad[0]].  
     174; left side ... if there is x2 values smaller that x[bad[0]]. 
    175175; we force ypn_1 = 0.0d 
    176176    IF xinx2[0] NE -1 THEN BEGIN 
    … …  
    178178        IF xinx2[0] NE 0 THEN stop 
    179179        y2[0] = y[0] 
    180       ENDIF ELSE BEGIN  
     180      ENDIF ELSE BEGIN 
    181181        y2[0:xinx2[0]] $ 
    182182          = spl_incr(x[0:bad[0]], y[0:bad[0]], x2[0:xinx2[0]] $ 
    183183                     , yp0 = yp0, ypn_1 = 0.0d) 
    184       ENDELSE  
    185     ENDIF  
     184      ENDELSE 
     185    ENDIF 
    186186; flat section 
    187187    IF xinx2_1[0] NE -1 THEN $ 
    … …  
    206206        ENDFOR 
    207207      ENDIF 
    208 ; right side ... if there is x2 values larger that x[bad[cntbad-1]+1].  
     208; right side ... if there is x2 values larger that x[bad[cntbad-1]+1]. 
    209209; we force yp0 = 0.0d 
    210210      IF xinx2_1[cntbad-1] NE nx2-1 THEN $ 
    … …  
    237237; 
    238238; we define the new values of the keyword ypn_1: 
    239 ; if the first derivative of the last value of x is negative  
     239; if the first derivative of the last value of x is negative 
    240240; we define the new values of the keyword ypn_1 to 0.0d0 
    241     IF bad[cntbad-1] EQ nx-1 THEN BEGIN  
     241    IF bad[cntbad-1] EQ nx-1 THEN BEGIN 
    242242      ypn_1new = 0.0d 
    243243; we remove this case from the list 
    … …  
    248248; 
    249249; we define the new values of the keyword yp0: 
    250 ; if the first derivative of the first value of x is negative  
     250; if the first derivative of the first value of x is negative 
    251251; we define the new values of the keyword yp0 to 0.0 
    252252    IF bad[0] EQ 0 THEN BEGIN 
    … …  
    265265; else: there is still cases with negative derivative ... 
    266266; we will cut spl_incr in n spl_incr and specify yp0, ypn_1 
    267 ; for each of this n spl_incr   
     267; for each of this n spl_incr 
    268268    ENDIF ELSE BEGIN 
    269269; define xinx2: see help of value_locate 
    … …  
    273273      xinx2 = value_locate(x2, x[bad]) 
    274274      y2 = dblarr(nx2) 
    275 ; left side ... if there is x2 values smaller that x[bad[0]].  
     275; left side ... if there is x2 values smaller that x[bad[0]]. 
    276276; we force ypn_1 = 0.0d 
    277277      IF xinx2[0] NE -1 THEN $ 
    … …  
    280280                        , yp0 = yp0new, ypn_1 = 0.0d) 
    281281; middle pieces ... if cntbad gt 1 then we have to cut spl_incr in 
    282 ; more than 2 pieces -> we have middle pieces for which  
     282; more than 2 pieces -> we have middle pieces for which 
    283283; we force yp0 = 0.0d and ypn_1 = 0.0d 
    284284      IF cntbad GT 1 THEN BEGIN 
    … …  
    295295        ENDFOR 
    296296      ENDIF 
    297 ; right side ... if there is x2 values larger that x[bad[cntbad-1]].  
     297; right side ... if there is x2 values larger that x[bad[cntbad-1]]. 
    298298; we force yp0 = 0.0d 
    299299      IF xinx2[cntbad-1] NE nx2-1 THEN $ 
    … …  
    302302                        , x2[xinx2[cntbad-1]+1:nx2-1] $ 
    303303                        , yp0 = 0.0d, ypn_1 = ypn_1new) 
    304     ENDELSE  
     304    ENDELSE 
    305305; we return the checked and corrected value of yfrst 
    306306;       FOR i = 0, nx-1 DO BEGIN 
    307307;         same = where(abs(x2- x[i]) LT 1.e-10, cnt) 
    308 ; ;        IF cnt NE 0 THEN y2[same] = y[i]  
     308; ;        IF cnt NE 0 THEN y2[same] = y[i] 
    309309;       ENDFOR 
    310310    RETURN, y2 
    … …  
    313313; we can be in this part of the code only if: 
    314314;  (1) spl_incr is called by itself 
    315 ;  (2) none are the first derivative in x are negative (because they have been  
    316 ;      checked and corrected by the previous call to spl_incr, see above)   
     315;  (2) none are the first derivative in x are negative (because they have been 
     316;      checked and corrected by the previous call to spl_incr, see above) 
    317317;--------------------------------- 
    318318; third check: we have to make sure that the first derivative cannot 
    … …  
    321321; 
    322322; first we compute the first derivative, next we correct the values 
    323 ; where we know that the first derivative can be negative.  
     323; where we know that the first derivative can be negative. 
    324324; 
    325325  y2 = spl_interp(x, y, yscd, x2, /double) 
    … …  
    330330; y''= 6a*X   + 2b 
    331331; if we take X = x[i+1]-x[i] then 
    332 ;    d = y[i]; c = y'[i]; b = 0.5 * y''[i],  
     332;    d = y[i]; c = y'[i]; b = 0.5 * y''[i], 
    333333;    a = 1/6 * (y''[i+1]-y''[i])/(x[i+1]-x[i]) 
    334 ;  
     334; 
    335335; y'[i] and y'[i+1] are positive so y' can be negative 
    336 ; between x[i] and x[i+1] only if  
     336; between x[i] and x[i+1] only if 
    337337;   1) a > 0 
    338338;            ==> y''[i+1] > y''[i] 
    339 ;   2) y' reach its minimum value between  x[i] and x[i+1]  
    340 ;      -> 0 < - b/(3a) < x[i+1]-x[i]  
     339;   2) y' reach its minimum value between  x[i] and x[i+1] 
     340;      -> 0 < - b/(3a) < x[i+1]-x[i] 
    341341;            ==> y''[i+1] > 0 > y''[i] 
    342342; 
    … …  
    412412; in those cases, the first derivative has 2 zero between 
    413413; x[bad[ib]] and x[bad[ib]+1]. We look for the minimum value of the 
    414 ; first derivative that correspond to the inflection point of y               
     414; first derivative that correspond to the inflection point of y 
    415415              xinfl = -bbb[ib]/(3.0d*aaa[ib]) 
    416416; we compute the y value for xinfl 
    417417              yinfl = aaa[ib]*xinfl*xinfl*xinfl + bbb[ib]*xinfl*xinfl $ 
    418418                + ccc[ib]*xinfl + ddd[ib] 
    419 ;                 
     419; 
    420420              CASE 1 OF 
    421421; if y[xinfl] smaller than y[bad[ib]] then we conserve y2 until 
    … …  
    450450                                     , yifrst[bad[ib]+1] $ 
    451451                                     , x2[xinx2_3+1:xinx2_2[ib]]) 
    452                   ENDIF                 
     452                  ENDIF 
    453453                END 
    454454; if y[xinfl] bigger than y[bad[ib]+1] then we conserve y2 from 
    … …  
    480480                                    , yifrst[bad[ib]] $ 
    481481                                    , x2[xinx2_1[ib]+1:xinx2_3]) 
    482                   ENDIF                 
     482                  ENDIF 
    483483                END 
    484484                ELSE:BEGIN 
    … …  
    496496                                    , yifrst[bad[ib]] $ 
    497497                                    , x2[xinx2_1[ib]+1:xinx2_3]) 
    498                      
    499                   ENDIF                 
     498 
     499                  ENDIF 
    500500                  IF xinx2_2[ib] GE xinx2_3+1 THEN BEGIN 
    501501                    y2[xinx2_3+1:xinx2_2[ib]] $ 
    … …  
    504504                                     , yifrst[bad[ib]+1] $ 
    505505                                     , x2[xinx2_3+1:xinx2_2[ib]]) 
    506                   ENDIF                 
     506                  ENDIF 
    507507                END 
    508508              ENDCASE 
    509509 
    510             END  
     510            END 
    511511          ENDCASE 
    512512        ENDIF 

  • trunk/SRC/Interpolation/spl_keep_mean.pro

    r121 r125  
    1414; data equa to the original values) 
    1515; 
    16 ;    @param x {in}{required}  An n-element (at least 2) input vector that specifies the 
    17 ;    tabulate points in a strict ascending order. 
     16; @param x {in}{required} 
     17; An n-element (at least 2) input vector that specifies the tabulate points in 
     18; a strict ascending order. 
    1819; 
    19 ;    @param yin {in}{required}  an array with one element less than x. y[i] represents the 
    20 ;    mean value between x[i] and x[i+1]. if /GE0 is activated, y must 
    21 ;    have positive values. 
     20; @param yin {in}{required} 
     21; an array with one element less than x. y[i] represents the 
     22; mean value between x[i] and x[i+1]. if /GE0 is activated, y must 
     23; have positive values. 
    2224; 
    23 ;    @param x2 {in}{required}  The input values for which the interpolated values are 
    24 ;    desired. Its values must be strictly monotonically increasing.  
     25; @param x2 {in}{required} 
     26; The input values for which the interpolated values are desired. 
     27; Its values must be strictly monotonically increasing. 
    2528; 
     29; @keyword GE0 
     30; to force that y2 is always GE than 0. In that case, y must also be GE than 0. 
    2631; 
    27 ; @keyword    /GE0 to force that y2 is always GE than 0. In that case, y must 
    28 ;    also be GE than 0. 
     32; @keyword YP0 
     33; The first derivative of the interpolating function at the 
     34; point X0. If YP0 is omitted, the second derivative at the 
     35; boundary is set to zero, resulting in a "natural spline." 
    2936; 
    30 ; @keyword    YP0 The first derivative of the interpolating function at the 
    31 ;    point X0. If YP0 is omitted, the second derivative at the 
    32 ;    boundary is set to zero, resulting in a "natural spline." 
     37; @keyword YPN_1 
     38; The first derivative of the interpolating function at the 
     39; point Xn-1. If YPN_1 is omitted, the second derivative at the 
     40; boundary is set to zero, resulting in a "natural spline." 
    3341; 
    34 ; @keyword    YPN_1 The first derivative of the interpolating function at the 
    35 ;    point Xn-1. If YPN_1 is omitted, the second derivative at the 
    36 ;    boundary is set to zero, resulting in a "natural spline."  
    37 ; 
    38 ; @returns  
     42; @returns 
    3943; 
    4044;    y2: the meean value between two consecutive values of x2. This 
    … …  
    4347; @restrictions 
    4448;   It might be possible that y2 has very small negative values 
    45 ;   (amplitude smaller than 1.e-6)...  
     49;   (amplitude smaller than 1.e-6)... 
    4650; 
    4751; 
    48 ; @examples  
     52; @examples 
    4953; 
    5054;    12 monthly values of precipitations into daily values: 
    … …  
    8993  nx2 = n_elements(x2) 
    9094; x must have at least 2 elements 
    91   IF nx LT 2 THEN stop  
     95  IF nx LT 2 THEN stop 
    9296; x2 must have at least 2 elements 
    93   IF nx2 LT 2 THEN stop  
     97  IF nx2 LT 2 THEN stop 
    9498; x be monotonically increasing 
    95   IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop  
     99  IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop 
    96100; x2 be monotonically increasing 
    97   IF min(x2[1:nx2-1]-x2[0:nx2-2])  LE 0 THEN stop  
     101  IF min(x2[1:nx2-1]-x2[0:nx2-2])  LE 0 THEN stop 
    98102; 
    99103;--------------------------------- 
    … …  
    130134;  yfrst = 0.0d > temporary(yfrst) 
    131135  RETURN, yfrst 
    132    
     136 
    133137;------------------------------------------------------------------ 
    134138;------------------------------------------------------------------ 

  • trunk/SRC/Interpolation/square2quadrilateral.pro

    r121 r125  
    11;+ 
    22; 
    3 ; @file_comments warm (or map) a unit square onto an arbitrary quadrilateral 
     3; @file_comments  
     4; warm (or map) a unit square onto an arbitrary quadrilateral 
    45; according to the 4-point correspondences: 
    56;       (0,0) -> (x0,y0) 
    … …  
    1314; @categories image, grid manipulation 
    1415; 
    15 ;  
    16 ;     @param x0in {in}{required}  the coordinates of the quadrilateral 
    17 ;     (see above for correspondance with the unit square). Can be 
    18 ;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
    19 ;     given in the anticlockwise order. 
    20 ;     @param y0in {in}{required}  the coordinates of the quadrilateral 
    21 ;     (see above for correspondance with the unit square). Can be 
    22 ;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
    23 ;     given in the anticlockwise order. 
    24 ;     @param x1in {in}{required}  the coordinates of the quadrilateral 
    25 ;     (see above for correspondance with the unit square). Can be 
    26 ;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
    27 ;     given in the anticlockwise order. 
    28 ;     @param y1in {in}{required}  the coordinates of the quadrilateral 
    29 ;     (see above for correspondance with the unit square). Can be 
    30 ;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
    31 ;     given in the anticlockwise order. 
    32 ;     @param x2in {in}{required}  the coordinates of the quadrilateral 
    33 ;     (see above for correspondance with the unit square). Can be 
    34 ;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
    35 ;     given in the anticlockwise order. 
    36 ;     @param y2in {in}{required}  the coordinates of the quadrilateral 
    37 ;     (see above for correspondance with the unit square). Can be 
    38 ;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
    39 ;     given in the anticlockwise order. 
    40 ;     @param x3in {in}{required}  the coordinates of the quadrilateral 
    41 ;     (see above for correspondance with the unit square). Can be 
    42 ;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
    43 ;     given in the anticlockwise order. 
    44 ;     @param y3in {in}{required}  the coordinates of the quadrilateral 
    45 ;     (see above for correspondance with the unit square). Can be 
    46 ;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
    47 ;     given in the anticlockwise order. 
    48 ; 
    49 ;     @param xxin {in}{optional} the coordinates of the point(s) for which we want to do the 
     16; @param x0in {in}{required}  the coordinates of the quadrilateral 
     17;     (see above for correspondance with the unit square). Can be 
     18;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
     19;     given in the anticlockwise order. 
     20; @param y0in {in}{required}  the coordinates of the quadrilateral 
     21;     (see above for correspondance with the unit square). Can be 
     22;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
     23;     given in the anticlockwise order. 
     24; @param x1in {in}{required}  the coordinates of the quadrilateral 
     25;     (see above for correspondance with the unit square). Can be 
     26;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
     27;     given in the anticlockwise order. 
     28; @param y1in {in}{required}  the coordinates of the quadrilateral 
     29;     (see above for correspondance with the unit square). Can be 
     30;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
     31;     given in the anticlockwise order. 
     32; @param x2in {in}{required}  the coordinates of the quadrilateral 
     33;     (see above for correspondance with the unit square). Can be 
     34;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
     35;     given in the anticlockwise order. 
     36; @param y2in {in}{required}  the coordinates of the quadrilateral 
     37;     (see above for correspondance with the unit square). Can be 
     38;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
     39;     given in the anticlockwise order. 
     40; @param x3in {in}{required}  the coordinates of the quadrilateral 
     41;     (see above for correspondance with the unit square). Can be 
     42;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
     43;     given in the anticlockwise order. 
     44; @param y3in {in}{required}  the coordinates of the quadrilateral 
     45;     (see above for correspondance with the unit square). Can be 
     46;     scalar or array. (x0,y0), (x1,y1), (x2,y2) and (x3,y3) are 
     47;     given in the anticlockwise order. 
     48; 
     49; @param xxin {in}{optional} the coordinates of the point(s) for which we want to do the 
    5050;     mapping. Can be scalar or array. 
    51 ;     @param yyin {in}{optional} the coordinates of the point(s) for which we want to do the 
     51; @param yyin {in}{optional} the coordinates of the point(s) for which we want to do the 
    5252;     mapping. Can be scalar or array. 
    5353; 
    … …  
    5555; 
    5656;     (2,n) array: the new coodinates (xout, yout) of the (xin,yin) 
    57 ;     point(s) after mapping.  
     57;     point(s) after mapping. 
    5858;     If xin is a scalar, then n is equal to the number of elements of 
    5959;     x0. If xin is an array , then n is equal to the number of 
    6060;     elements of xin. 
    6161;     If xin and yin are omited, square2quadrilateral returns the 
    62 ;     matrix A which is used for the inverse transformation.  
     62;     matrix A which is used for the inverse transformation. 
    6363; 
    6464; 
    … …  
    6666; twisted) is not work. This has to be tested. 
    6767; 
    68 ; @examples  
     68; @examples 
    6969; 
    7070; IDL> splot,[0,5],[0,3],/nodata,xstyle=1,ystyle=1 
    … …  
    8181;      IEEE Computer Society Press, Los Alamitos, California 
    8282;      Chapter 3, see p 52-56 
    83 ;       
     83; 
    8484; 
    8585; @version $Id$ 
    … …  
    127127; 
    128128  IF keyword_set(double) THEN a = dblarr(8, n_elements(x0)) $ 
    129   ELSE a = fltarr(8, n_elements(x0))  
     129  ELSE a = fltarr(8, n_elements(x0)) 
    130130; 
    131131  delx3 = x0-x1+x2-x3 
    … …  
    161161    yy2 = y2[projectivemap] 
    162162    yy3 = y3[projectivemap] 
    163 ;     
     163; 
    164164    delx1 = xx1-xx2 
    165165    dely1 = yy1-yy2 
    … …  
    186186    a[7, projectivemap] = a23 
    187187  ENDIF 
    188 ;     
     188; 
    189189  IF NOT arg_present(xxin) THEN return, a 
    190190; 
    191191  IF n_elements(xin) EQ 1 THEN BEGIN 
    192     xin = replicate(xin, n_elements(x0))  
    193     yin = replicate(yin, n_elements(x0))  
     192    xin = replicate(xin, n_elements(x0)) 
     193    yin = replicate(yin, n_elements(x0)) 
    194194  ENDIF 
    195195; 
Note: See TracChangeset for help on using the changeset viewer.