FUNCTION gridgen, nn, yy, zz, $ unit = unit , $ merge = merge , $ origin = origin , $ edge = edge , $ open = open , $ one = one , $ range = range , $ multid = multid , $ double = double @compile_opt.pro ; On error, return to caller InitVar, multid, /key ; This is a bit sick, but I have done worse InitVar, merge, /key IF merge THEN BEGIN CASE n_params() OF 1: r = nn 2: BEGIN n1 = n_elements(nn) n2 = n_elements(yy) r = make_array(dim=[2,n1*n2], type=(IsType(nn) > IsType(yy)), /nozero) r[0,*] = SuperArray(nn,n2,/trail) r[1,*] = SuperArray(yy,n1,/lead ) IF multid THEN r = reform(r,2,n1,n2,/overwrite) ELSE r = reform(r,/overwrite) END 3: BEGIN n1 = n_elements(nn) n2 = n_elements(yy) n3 = n_elements(zz) r = make_array(dim=[3,n1*n2*n3], type=(IsType(nn) > IsType(yy) > IsType(zz)), /nozero) r[0,*] = SuperArray(nn,n2*n3,/trail) r[1,*] = SuperArray(SuperArray(yy,n1,/lead),n3,/trail) r[2,*] = SuperArray(zz,n1*n2,/lead) IF multid THEN r = reform(r,3,n1,n2,n3,/overwrite) ELSE r = reform(r,/overwrite) END ELSE: message, '/merge can only handle 1-,2- and 3-D grids' ENDCASE RETURN, r ENDIF IF IsType(range ,/string) THEN range = TimeSet(range ) IF IsType(origin,/string) THEN origin = TimeSet(origin) time_grid = IsTime(range) OR IsTime(origin) InitVar, nn, 0 ; nn may be absent n = long(nn) ndim = n_elements(n) multid AND= ndim GT 1 has_origin = IsType(origin , /defined) has_edge = IsType(edge , /defined) has_open = IsType(open , /defined) has_one = IsType(one , /defined) has_range = IsType(range, /defined) IF has_range THEN BEGIN range_ = range IF n[0] EQ 0 THEN BEGIN ; No # elements specified CASE time_grid OF 0: BEGIN n = n_elements(range_) CASE 1 OF n EQ 1 : range_ = [0,range_] n/2*2 EQ n: range_ = reform(range_,2,n/2,/overwrite) ELSE : message, 'illegal range spec' ENDCASE InitVar, unit, 1 n = round(reform((range_[1,*]-range_[0,*]))/unit)+1 ndim = n_elements(n) multid AND= ndim GT 1 END 1: BEGIN InitVar, unit, TimeUnit(/day) n = round(TimeOp(/subtract,range_[1],range_[0],unit))+1 range_ = range_[0] END ENDCASE ENDIF CASE n_elements(range_) OF 1 : BEGIN CASE time_grid OF 0: range_ = [0,range_] 1: BEGIN InitVar, unit, TimeUnit(/day) range_ = [range_,TimeOp(/add,range_,TimeSet(n-1,unit,/diff))] END ENDCASE END ndim: range_ = [intarr(1,ndim), reform(range_,1,ndim)] ELSE: range_ = reform(range_,2,ndim,/overwrite) ENDCASE ENDIF n_origin = n_elements(origin) -1 n_edge = n_elements(edge ) -1 n_open = n_elements(open ) -1 n_one = n_elements(one ) -1 n_range = n_elements(range_)/2-1 ; The original gridgen went upto ndim=4 using the following lines: ; r[0,*] = reform( reform( gridgen1d(n[0],origin=origin[0],edge=edge[0],open=open[0],one=one[0],range=range_[0:1],double=double) # replicate(1L,n[1]), n01 ) # replicate(1L,n[2]), n012 ) # replicate(1L,n[3]) ; r[1,*] = reform( reform( replicate(1L,n[0]) # gridgen1d(n[1],origin=origin[1],edge=edge[1],open=open[1],one=one[1],range=range_[2:3],double=double), n01 ) # replicate(1L,n[2]), n012 ) # replicate(1L,n[3]) ; r[2,*] = reform( reform( replicate(1L,n[0]) # replicate(1L,n[1]), n01 ) # gridgen1d(n[2],origin=origin[2],edge=edge[2],open=open[2],one=one[2],range=range_[4:5],double=double), n012 ) # replicate(1L,n[3]) ; r[3,*] = reform( reform( replicate(1L,n[0]) # replicate(1L,n[1]), n01 ) # replicate(1L,n[2]), n012 ) # gridgen1d(n[3],origin=origin[3],edge=edge[3],open=open[3],one=one[3],range=range_[6:7],double=double) ; Calculate number of elements in grid ; Remember that gridgen1d adds an element if keyword edge is set. np = lonarr(ndim) IF has_edge THEN np[0] = n[0]+edge[0] ELSE np[0] = n[0] FOR i=1,ndim-1 DO IF has_edge THEN np[i] = np[i-1]*(n[i]+edge[i < n_edge]) ELSE np[i] = np[i-1]*n[i] FOR i=0,ndim-1 DO BEGIN ; Loop over all ndim components FOR j=0,i-1 DO BEGIN grid = replicate(1L,n[j]) IF j EQ 0 THEN tmp = grid ELSE tmp = reform( tmp # grid, np[j] ) ENDFOR IF has_origin THEN origin_i = origin[ i < n_origin] IF has_edge THEN edge_i = edge [ i < n_edge ] IF has_open THEN open_i = open [ i < n_open ] IF has_one THEN one_i = one [ i < n_one ] IF has_range THEN range_i = range_[*,i < n_range ] grid = gridgen1d(n[i],unit,origin=origin_i,edge=edge_i,open=open_i,one=one_i,range=range_i,double=double) CASE i NE 0 OF 0: BEGIN tmp = grid ; The first call to gridgen1d sets the type of the output array. CASE time_grid OF 0: r = make_array(type=size(grid,/type), dimension=[ndim,np[ndim-1]]) 1: r = replicate( !TheTime, [ndim,np[ndim-1]] ) ENDCASE ENDCASE 1: tmp = reform(tmp # grid, np[i]) ENDCASE FOR j=i+1,ndim-1 DO BEGIN grid = replicate(1L,n[j]) tmp = reform( tmp # grid, np[j] ) ENDFOR r[i,*] = tmp ENDFOR IF multid THEN r = reform(r, [ndim,n], /overwrite) ELSE r = reform(r, /overwrite) RETURN, r & END