Changeset 40 for trunk


Ignore:
Timestamp:
06/10/11 06:00:38 (12 months ago)
Author:
tabel
Message:

a few bug fixes

Location:
trunk
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • trunk/TOOLS/AMRslice.pro

    r39 r40  
    616616           step = 1 
    617617           exists = '' 
    618            while ((step lt 1000) and (exists[0] eq '')) do begin 
     618           while ((step lt 200) and (exists[0] eq '')) do begin 
    619619;              findS = sp[0:ns-2] + STRING(max([LONG(timenum)-step,0]), FORMAT='(i4.4)') 
    620620              if regular then $ 
    … …  
    672672           sp = strsplit(file_name, timenum, /extract, /regex) 
    673673           ns = N_elements(sp) 
    674            step = 1 
     674           step = 0 
    675675           exists = '' 
    676            while ((step lt 1000) and (exists[0] eq '')) do begin 
     676           while (step le astep) do begin 
     677           exists = '' 
     678           while ((step lt 200) and (exists[0] eq '') ) do begin 
    677679              if regular then $ 
    678680                 findS = sp[0:((ns-2) > 0)] + STRING(LONG(timenum)+step, FORMAT='(i4.4)') $ 
    … …  
    682684              protofname = strjoin(findS) 
    683685              exists = findfile(protofname) 
     686           endwhile 
    684687           endwhile 
    685688           if exists[0] ne '' then begin 
    … …  
    973976              if last_was eq 1 then j_event, 'DRAW_BUTTON'  
    974977              if last_was eq 2 then j_event, 'PROJECT' 
    975               if last_was eq 3 then j_event, 'DRAW_BUTTON' 
     978              if last_was eq 3 then j_event, 'PROJECT' 
    976979              if last_was eq 4 then j_event, 'PROJECTABS' 
    977980              if last_was eq 5 then j_event, 'PROJECTHII' 
    … …  
    988991              if last_was eq 1 then j_event, 'DRAW_BUTTON'  
    989992              if last_was eq 2 then j_event, 'PROJECT' 
    990               if last_was eq 3 then j_event, 'DRAW_BUTTON' 
     993              if last_was eq 3 then j_event, 'PROJECT' 
    991994              if last_was eq 4 then j_event, 'PROJECTABS' 
    992995              if last_was eq 5 then j_event, 'PROJECTHII' 
    … …  
    10601063           construct_histogram, data 
    10611064           last_was = 3 
     1065           draw_current_cm, draw_cm, wxysize, 0., 1., info_only=1 
    10621066        endif else if parti.myturn eq 1 then begin 
    10631067           data = render_particle_data() 
    … …  
    12811285 
    12821286        show_hierarchy_information, grid_info 
     1287        print, 'total number of particles:', total(grid_info.num_particle) 
    12831288    END 
    12841289     
    … …  
    17061711;        if useunits then unitfactor = get_unit(image.title, /ufactor) 
    17071712        show_image 
    1708         draw_current_cm, draw_cm, wxysize, image_stack[image_index].min_data, image_stack[image_index].max_data 
     1713;        draw_current_cm, draw_cm, wxysize, image_stack[image_index].min_data, image_stack[image_index].max_data 
    17091714    END 
    17101715 

  • trunk/TOOLS/OJ__define.pro

    r1 r40  
    10751075   endif else hdat = get_data(use_g, data_fields, viewdim,/return_as_list) 
    10761076    
    1077  
    10781077   n_df_requested = N_elements(data_fields) 
    10791078 
    10801079; how many data points in total? 
    1081    tot_cells = hdat[0].np_total 
     1080   tot_cells = N_elements(*hdat[0].data[0]) > 1 
    10821081; allocate some memory 
    10831082   if n_df_requested ge 1 then d_stuetz = DBLARR(N_df_requested,tot_cells) 
    … …  
    19061905    num_pattr = 8             ; # of Particle Attributes:  
    19071906                               ; Star formation = 11, None = 8 
    1908     max_num_of_grids = 150000L 
     1907    max_num_of_grids = 250000L 
    19091908    max_num_of_levels = 45 
    19101909 

  • trunk/TOOLS/construct_amira_export.pro

    r6 r40  
    4242  last_slash = STRPOS(parts[0], "/", /REVERSE_SEARCH)+1 
    4343  filename = STRMID(parts[0], last_slash) 
    44   return, [filename, parts[1]] 
     44  return, [filename, parts[0]] 
    4545 
    4646end 

  • trunk/TOOLS/construct_histogram.pro

    r36 r40  
    137137        mh = hist1d(x, $ 
    138138                        min=min1, max=max1, $ 
    139                         binsize=(max1-min1)/hi.nbins1, obin=obin1,DENSITY=dens) 
     139                        binsize=(max1-min1)/(hi.nbins1-1), obin=obin1,DENSITY=dens) 
    140140        if (hi.cumoption eq 1) then mh = total(mh, /cumulative) 
    141141        if (hi.cumoption eq 2) then mh = reverse(total(reverse(mh), /cumulative)) 
    … …  
    143143        if ((not useunits) and (hi.cumoption gt 0)) then mh = 1.D*mh/max(mh)  
    144144 
     145        obin1 = obin1- 0.5*(max1-min1)/N_elements(obin1) 
    145146        if shouldbe_logged(validfields[0]) then begin 
    146147           obin1=10.^obin1 
    … …  
    151152        if hi.max2 eq  1e30 then max2 = max(mh) else max2 = hi.max2 
    152153 
    153         print ,'show hist for one', nnotnone, obin1, mh 
     154;        print ,'show hist for one', nnotnone, obin1, mh 
    154155        plot, obin1, mh, xlog=shouldbe_logged(validfields[0]), $ 
    155156              ylog=hi.logoption, psym=10, $ 
    … …  
    177178           if hi.max2 eq  1e30 then max2 = max(mh) else max2 = hi.max2 
    178179            
    179            print ,'show hist for two', nnotnone, obin1, mh 
     180;           print ,'show hist for two', nnotnone, obin1, mh 
    180181           plot, obin1, (mh) , xlog=shouldbe_logged(validfields[0]), $ 
    181182                 ylog=hi.logoption, psym=10, $ 
    … …  
    267268  hi.xrange = [min(obin1),max(obin1)] 
    268269  hi.yrange  =[min(obin2),max(obin2)] 
    269    
     270 
    270271  x_stuetz = x 
    271272  if Nnotnone gt 1 then y_stuetz = y 
    272273  if Nnotnone gt 2 then z_stuetz = z 
     274 
     275  if shouldbe_logged(validfields[0]) then x_stuetz = 10.^x 
     276  if (Nnotnone gt 1) then if (notnone[1] eq 1) then  if shouldbe_logged(validfields[1]) then y_stuetz = 10.^y 
     277  if Nnotnone gt 2 then  if shouldbe_logged(validfields[2]) then z_stuetz = 10.^z 
    273278   
    274279  return 

  • trunk/TOOLS/construct_interpolated_slice_data.pro

    r38 r40  
    234234 
    235235  n_stuetz = 0L 
    236   FOR i=0L,N_elements(a)-1 DO BEGIN 
     236  FOR i=0L, N_elements(a)-1 DO BEGIN 
     237;  FOR i=0L, 0 DO BEGIN 
    237238; check this 
    238239     i_s   = ii_s[i] 

  • trunk/TOOLS/draw_current_cm.pro

    r34 r40  
    124124END ;---------------------------------------------------------------------------------------- 
    125125 
    126 PRO draw_current_cm, draw_widget_id, x_size, min_v, max_v, return_strings=rstrings 
     126PRO draw_current_cm, draw_widget_id, x_size, min_v, max_v, return_strings=rstrings, info_only=info_only 
    127127@common_blocks.inc 
    128128   image = image_stack[image_index] 
    … …  
    199199   im_arr[x_off:(x_off +width-1), y_off:(y_off+height-1)] = $ 
    200200      BYTSCL(CONGRID(INDGEN(ncolors-2),width,/MINUS_ONE,/INTERP,/CENTER) # REPLICATE(1, height)) + 2 
    201    if last_was ne 6 then begin  ; draw color bar 
     201   if ((last_was ne 6) and NOT KEYWORD_SET(info_only)) then begin  ; draw color bar 
    202202      TVSCL, im_arr 
    203203;draw bounding box 
    … …  
    247247    
    248248   if last_was ne 0 then titlestring=logstring+image.title+unitstring else titlestring = 'AMR Level'  
    249    xyouts, 10, 55, titlestring, /DEVICE, CHARSIZE=2., CHARTHICK=1, COLOR=ncolors-1 
     249   if (NOT KEYWORD_SET(info_only)) then $ 
     250      xyouts, 10, 55, titlestring, /DEVICE, CHARSIZE=2., CHARTHICK=1, COLOR=ncolors-1 
    250251;device, set_font='9x15bold' 
    251252   if use_nice_font then device, set_font='-*-helvetica-*-r-*-*-*-160-*-*-*-*-*-*' 

  • trunk/TOOLS/draw_grids.pro

    r27 r40  
    3434      delta_dist   = delta_distance(*,i) 
    3535      data_points  = index_range(*,i)   
    36       slice_ori_l                = DBLarr(3) 
     36      slice_ori_l                = DBLARR(3) 
    3737      slice_ori_l(const_sub)     = slice_ori(const_sub) 
    3838      slice_ori_l(other_subs)    = [slice_coord(0),slice_coord(1)] 
    … …  
    4040      slice_ori_r(other_subs)    = [slice_coord(2),slice_coord(3)] 
    4141 
    42       i_s = ROUND((slice_ori_l(const_sub) - $ 
     42      i_s = FLOOR((slice_ori_l(const_sub) - $ 
    4343                   u_g(i).Left_edge(const_sub))/delta_dist(const_sub) >0 $ 
    4444        < data_points(const_sub)-1) 
    45       i_l = CEIL( (slice_ori_l(other_subs) - $ 
     45      i_l = FIX( (slice_ori_l(other_subs) - $ 
    4646                   u_g(i).Left_edge(other_subs)) $ 
    4747                    /delta_dist(other_subs) ) > 0 
    48       i_r   = (FLOOR((slice_ori_r(other_subs) - $ 
     48      i_r   = (FIX((slice_ori_r(other_subs) - $ 
    4949                  u_g(i).Left_edge(other_subs))/ $ 
    50                     delta_dist(other_subs) - 1) < $ 
     50                    delta_dist(other_subs) ) < $ 
    5151        (data_points(other_subs)-1) ) > 0 
    5252 

  • trunk/TOOLS/export_particles_for_ralf.pro

    r1 r40  
    3939  nfields = N_elements(names) 
    4040  printf, lun, '#point format(ASCII) V00' 
    41   for i=0, 3 do begin 
     41  for i=3, 3 do begin 
    4242     thisi = i 
    4343     if i eq 3 then thisi = parti.index  ; make sure the 4th colums is what is selected in GUI 
    44      printf, lun, '#datavar'+strcompress(i-3)+' '+strcompress(names[thisi], /remove_all) 
     44     printf, lun, '#fieldnames: '+strcompress(names[thisi], /remove_all) 
    4545  endfor 
    4646 
    … …  
    8484           endfor 
    8585        endfor 
    86         printf, lun, '' 
     86;        printf, lun, '' 
    8787     endif 
    8888     close, lun 

  • trunk/TOOLS/find_max_value.pro

    r17 r40  
    5858          index(2)   =  full_index/ $ 
    5959            (data_points(0) * data_points(1)) mod data_points(2) 
    60           xyz_vector = DOUBLE(grid_info[i].Left_edge) + (index)*delta_dist 
     60          xyz_vector = DOUBLE(grid_info[i].Left_edge) + (index+0.5)*delta_dist 
    6161          print, 'xyz_vec, max, maxcheck:', xyz_vector, max_value,$ 
    6262            (*a[i].data[0])[index(0),index(1),index(2)] 

  • trunk/TOOLS/get_data.pro

    r39 r40  
    99pro pack_as_list, use_g, a 
    1010; pack all data into 1D arrays  
    11    zero = -1e30 
     11   zero = 0. 
    1212   zero_solution_under_subgrid, use_g, a, ZERO=zero 
    1313 
    … …  
    359359           ra = gi[j].end_index-gi[j].start_index + 1 
    360360           tmp = fltarr(ra[0], ra[1], ra[2]) 
    361            if verbose then print, 'debug:', ra, N_elements(tmp) 
     361;           if verbose then print, 'debug:', ra, N_elements(tmp) 
    362362           tmp[*] = 1. 
    363363           a[j].li = [0, 0, 0] 
    … …  
    402402                  ri[k] = ((min([gi[j].parent_e_index[k]+1, a[pin].ri[k]])-a[pin].li[k]) $ 
    403403                           < (a[pin].dims[k])-1) > (li[k] + 1) 
    404                print, 'li,ri:', li, ri 
    405                print, 'd', d 
    406                print, 'a[pin].dims:',a[pin].dims 
    407                print, 'a[pin].li/ri:',a[pin].li, a[pin].ri 
     404;               print, 'li,ri:', li, ri 
     405;               print, 'd', d 
     406;               print, 'a[pin].dims:',a[pin].dims 
     407;               print, 'a[pin].li/ri:',a[pin].li, a[pin].ri 
    408408               if min(d) gt 2  then begin 
    409409; should be doing something like the following ... however it does not 
    … …  
    518518         17: *a[j].data[i] = fn(a[j], field_names[0]) * fn(a[j], "Density") 
    519519         18: *a[j].data[i] = sqrt(!DPI*3./(32.*fn(a[j], "Density"))) 
    520          19: *a[j].data[i] = sqrt((center[0]-fn(a[j], "x"))^2+(center[1]-fn(a[j], "y"))^2 $  ; radius 
    521                                   +(center[2]-fn(a[j], "z"))^2) 
     520         19: *a[j].data[i] = sqrt((center[0]-fn(a[j], "x"))^2+(center[1]-fn(a[j], "y"))^2 $ ; radius 
     521                                 +(center[2]-fn(a[j], "z"))^2) 
    522522         20: *a[j].data[i] = sqrt(rcd[0]*(center[0]-fn(a[j], "x"))^2 + $ 
    523523                                  rcd[1]*(center[1]-fn(a[j], "y"))^2 + $  
    … …  
    539539         28: begin 
    540540            uvel = get_unit("Velocity", /force_use_units) 
    541             *a[j].data[i] = 1.1725e-1 * sqrt(fn(a[j], "Temperature")) / uvel 
     541            *a[j].data[i] = 1.1725e-1 * sqrt(fn(a[j], "Temperature")) * sqrt(gamma/1.6667) / uvel 
    542542         end 
    543543         29: begin 
    544544            umass = get_unit("mass", /force_use_units) 
    545             *a[j].data[i] = 19.0315 * (fn(a[j], "Temperature"))^1.5 / $ 
     545            *a[j].data[i] = 19.0315 * ((gamma/1.6667) * fn(a[j], "Temperature"))^1.5 / $ 
    546546                            sqrt(fn(a[j], "Density")) / umass 
    547547         end 
    … …  
    557557                              DOUBLE(fn(a[j], "x-velocity"))^2 + $ 
    558558                              DOUBLE(fn(a[j], "y-velocity"))^2 + $ 
    559                               DOUBLE(fn(a[j], "z-velocity")^2 )));*fn(a[j], "Density") 
     559                              DOUBLE(fn(a[j], "z-velocity")^2 )))*fn(a[j], "Density") 
    560560 
    561561         35: *a[j].data[i] =  ALOG(1.D*(0.58D/0.6D*fn(a[j], "Temperature")) /((1.+z)^3*fn(a[j], "Density")*6.9416d10)^(gamma-1.)) ; entropy in funny units ... assuming was mu=0.6 in enzo sim 
    562562 
    563          36: *a[j].data[i] = div_v(fn(a[j], "x-velocity"),fn(a[j], "y-velocity"),fn(a[j], "z-velocity"), $ 
    564                                    fn(a[j], "dx"),fn(a[j], "dy"),fn(a[j], "dz")) 
     563         36: *a[j].data[i] = (div_v(fn(a[j], "x-velocity"),fn(a[j], "y-velocity"),fn(a[j], "z-velocity"), $ 
     564                                   fn(a[j], "dx"),fn(a[j], "dy"),fn(a[j], "dz"))) 
    565565           
    566566         37: *a[j].data[i] = rho_slope(fn(a[j], "Density"),fn(a[j], "x"), $ 

  • trunk/TOOLS/get_units.pro

    r39 r40  
    129129        unitfactor = velocityunit/1e5/(1.+InitialRedshift) 
    130130        unitlabel = 'km/s' 
     131        if (data_format eq 7) and (label eq "Alfven speed") then begin 
     132            print, 'flash format! ' 
     133            unitfactor /= sqrt(4.D*!DPi) 
     134        endif 
    131135     end 
    132136     (tlabel eq 'Energy'): Begin 
    … …  
    155159;        unitfactor = sqrt(densityunit*(velocityunit)^2*4.*!pi)/1.e-6 
    156160;        unitlabel = 'micro Gauss' 
    157         unitfactor =  (1.D + InitialRedshift)^0.5 * sqrt(densityunit* $ 
    158                           (velocityunit)^2*4.*!pi) * $ 
    159                      ((1.D + InitialRedshift)/(1.+zz))^(-1.5) 
     161 
     162        unitfactor =  sqrt(densityunit*(1.+zz)^3 *(velocityunit)^2 * 4.*!pi)  
     163 
    160164        unitlabel = 'Gauss' 
     165 
     166        if (data_format eq 7) then begin ;  flash differs by factor 4 Pi in magnetic units 
     167            print, 'flash format! ' 
     168            unitfactor /= sqrt(4.D*!DPi) 
     169        endif 
     170 
    161171     end 
    162172     (tlabel eq 'mass'): Begin 
    … …  
    170180     end 
    171181     (tlabel eq '|v_particles|^10'): begin 
    172         unitfactor = (velocityunit/1e5/(1.+InitialRedshift))^10 
     182        unitfactor = (velocityunit/1e5)^10 
    173183        unitlabel = textoidl('(km s^{-1})^{10}') 
    174184     end 
    175185     (tlabel eq 'TT12'): begin 
    176         unitfactor = (velocityunit/1e5/(1.+InitialRedshift))^4 * (densityunit*lengthunit^3/1.989d33)^2 
     186        unitfactor = (velocityunit/1e5)^4 * (densityunit*lengthunit^3/1.989d33)^2 
    177187        unitlabel = textoidl('(km s^{-1})^{4} ')+textoidl('M')+'!D!9n!X!N'+textoidl('^2') 
    178188     end 
    … …  
    212222        unitfactor = magunitfactor/lengthunit 
    213223        unitlabel  = 'Gauss/cm' 
     224        if (data_format eq 7) then begin ;  flash differs by factor 4 Pi in magnetic units 
     225            print, 'flash format! ' 
     226            unitfactor /= sqrt(4.D*!DPi) 
     227        endif 
     228 
    214229     end 
    215230 

  • trunk/TOOLS/grid_in_cube.pro

    r39 r40  
    4141 
    4242  print, 'there are ', total(help),' grids of sufficient resolution at least partially in cube' 
    43   if total(help) gt 1 then grid_in = new_i(where(help gt 0)) 
     43  if total(help) ge 1 then grid_in = new_i(where(help gt 0)) 
    4444  return, grid_in 
    4545END 

  • trunk/TOOLS/hdf5_read_routines.pro

    r34 r40  
    251251 
    252252           if (where(field eq special))[0] ne -1 then begin ; this is a special field 
    253                if count_defined then ra = ((count[cnt,*]) > 1 ) else ra=cg.dim 
     253               if count_defined then ra = ((count[cnt,*]) > 1 ) else ra=cg.dim > 1 
    254254               in = INDGEN(product(ra), /l64) 
    255255               data = dblarr(ra[0],ra[1],ra[2]) 
    256256;print,ra, count[cnt,*] 
     257;               print, 'ra',ra ,'dim',cg.dim 
    257258               indi = array_indices(ra,in, /dimensions) 
    258259               case field of 
    … …  
    263264                   'dy': data[*] = dd[1] 
    264265                   'dz': data[*] = dd[2] 
    265                    'volume': data[*] = product(dd[where(dims gt 0)]) 
     266                   'volume': data[*] = product(dd[where(cg.dim gt 1)]) 
    266267                   ELSE: data[*] = 1. 
    267268               endcase 
    268  
    269269;debs= size(data) 
    270270;if debs[0] ne 3 then print, debs, ra 
    … …  
    274274                 hdfgid = H5G_OPEN(fid, gident) 
    275275                 ngroup = H5G_GET_NUM_OBJS(hdfgid) 
     276                 H5G_CLOSE, hdfgid 
    276277                 kl = [''] 
    277278                 for gm=0,ngroup-1 do begin 
    … …  
    312313                   H5D_CLOSE, dsetid 
    313314               endif else print, 'warning requested field not found in ', gident, field 
    314            ENDELSE      
     315            ENDELSE 
    315316            
    316317           all[cnt].file_name   = gident 
    317318           all[cnt].np_total    = N_elements(data) 
    318319           s = size(data) 
    319            if (s[0] eq 2) then s[3] = 1 
     320           if (s[0] le 2) then s[3] = 1 
     321           if (s[0] eq 1) then s[2] = 1 
    320322           all[cnt].li[*]       = this_start[*] 
    321323           all[cnt].ri[*]       = (s[1:3]+this_start-1)[*] 

  • trunk/TOOLS/make_hierarchy_from_movie_format_file.pro

    r1 r40  
    5959 
    6060     ;; Number of grids is written at the end of the file 
     61     if index_name eq '' then break 
    6162     openr, 1, index_name 
    6263     filesize = (FSTAT(1)).size 
    … …  
    120121  endfor 
    121122 
     123  if N_elements(grid_info) eq 0 then return, '' 
    122124  grid_info = grid_info[sort(grid_info.time)] 
    123125 

  • trunk/TOOLS/read_grid_info.pro

    r29 r40  
    8888  temp_grid_info = info 
    8989 
    90   if strmatch(file_name, '*.hierarchy') then BEGIN 
     90  nlines = 0 
     91  if strmatch(file_name, '*.hierarchy') then $ 
     92     fstring = read_ascii_file_return_list_of_strings(file_name[0], /verbose, nlines=nlines)  
     93  if (nlines gt 0) then BEGIN 
    9194  ; hdf4 or packed hdf5 format 
    92      fstring = read_ascii_file_return_list_of_strings(file_name[0], /verbose, nlines=nlines)  
    9395     dim_lines = fstring[where((strmatch(fstring, 'GridDim*') eq 1), ngrids)] 
    9496     temp_grid_info  = Replicate(info, ngrids)  
    … …  
    114116     ind = LONARR(N_elements(hngnl)) 
    115117     equalpos = strpos(hngnl,"=") 
    116      for ii=0L, ngrids-1 do BEGIN 
     118     for ii=0L, N_elements(hngnl)-1 do BEGIN 
    117119        mi = LONG(STRMID(hngnl[ii], 14, strpos(hngnl[ii], ']')-14))-1 
    118120        mngnl = LONG(strmid(hngnl[ii], equalpos[ii]+1)) 
    … …  
    125127     ind = LONARR(ngrids) 
    126128     equalpos = strpos(hngtl,"=") 
    127      for ii=0L, ngrids-1 do $ 
     129     for ii=0L, N_elements(hngtl)-1 do $ 
    128130        ngtl(LONG(STRMID(hngtl[ii], 14, strpos(hngtl[ii], ']')-14))-1) = long(strmid(hngtl[ii], equalpos[ii]+1)) 
    129131      

  • trunk/TOOLS/read_parameters.pro

    r35 r40  
    9999;        mpc*ComovingBoxSize/HubbleConstantNow/(1.+redshift) 
    100100        lengthunit = mpc*ComovingBoxSize/HubbleConstantNow 
     101        velocityunit = lengthunit/timeunit/(1.+InitialRedshift) 
    101102     endif 
    102103      
    … …  
    150151     endif 
    151152 
    152      velocityunit = lengthunit/timeunit 
     153     if (InitialRedshift eq 0.) then velocityunit = lengthunit/timeunit 
    153154 
    154155end 

  • trunk/TOOLS/read_velocities.pro

    r1 r40  
    1010; magnetic field or velocities? 
    1111      if (vector_field_b GT 0) then check = STRPOS(list_str(i),'Bx') else $ 
    12                check = STRPOS(list_str(i),'elocity')  
     12               check = STRPOS(list_str(i),'-velocity')  
    1313      pcheck = STRPOS(list_str(i),'articl')  
    1414      if ((check ge 0) and (pcheck eq -1) and (found_at eq 0)) THEN found_at = i 

  • trunk/TOOLS/shouldbe_logged.pro

    r34 r40  
    99         (strpos(string, 'orticity') GE 0) OR $ 
    1010         (strpos(string, 'Entropy') GE 0) OR $ 
     11         (strpos(string, 'AccelerationField') GE 0) OR $ 
    1112         (strpos(string, 'PotentialField') GE 0) OR $ 
    1213         (strpos(string, 'Grav_Potential') GE 0) OR $ 

  • trunk/TOOLS/zero_solution_under_subgrid.pro

    r1 r40  
    1414                ri[j] = ((min([a_g[cc].parent_e_index[j], all[i].ri[j]])-all[i].li[j]) $ 
    1515                         < (all[i].dims[j]-1)) > li[j] 
    16 ;              print, ri-li+1 
     16;              print, 'zsus:', ri-li+1 
    1717              for j=0,all[i].N_Dfields-1 do begin 
    1818                 if all[i].np_total gt 0 then begin 
    1919                    ds = size(*all[i].data[j]) 
    20                     if ds[0] lt 2 then break 
    21                     if N_elements(*all[i].data[j]) eq product(all[i].dims > 1) then $ 
    22                        (*all[i].data[j])[li[0]:ri[0], li[1]:ri[1], li[2]:ri[2]] = zero else begin 
     20                    if ds[0] lt 1 then break 
     21                    if N_elements(*all[i].data[j]) eq product(all[i].dims > 1) then begin 
     22                       (*all[i].data[j])[li[0]:ri[0], li[1]:ri[1], li[2]:ri[2]] = zero  
     23;                       print, li, ri 
     24                    endif else begin 
    2325                       print, 'zero_solution_under_subgrid: somethings wrong', a_g[cc].num, i, jj, j, $ 
    24                               all[i].data_fields[j],N_elements(all[i].data[j]), product(all[i].dims),(all[i].dims) 
     26                              all[i].data_fields[j], 'Ne:', N_elements(all[i].data[j]), product(all[i].dims),(all[i].dims) 
     27;                       breaki 
    2528                    endelse 
    2629                 endif 

  • trunk/scripts/time_averaged_density_pdf.pro

    r30 r40  
    1818 
    1919  tapdf = fltarr(hi.nbins1)     ; field to store result 
     20  tclump = 0. 
    2021  totaltime = 0. 
    2122 
    … …  
    3233     wait, 1. 
    3334     select_current_grids       ; 
    34      check_cancel              ; allow user to type  c in command line to stop  
     35     if check_cancel() then return ; allow user to type  c in command line to stop  
    3536     if stopit then begin 
    3637        stopit = 0 
    … …  
    4344     tapdf += mh 
    4445 
     46     clump = mean(x*x)/mean(x)^2 
     47     tclump += clump 
     48     print, 'clumping factor:' , clump 
     49      
    4550     if verbose then print, 'memory use in Mb', memory()/1024./1e3 
    4651  endfor 
    … …  
    4853 
    4954  tapdf /= totaltime            ; compute the average 
     55  tclump /= totaltime           ; compute average clumping factor 
    5056 
    5157  mh = tapdf                    ;  store in this common block variable 
    … …  
    6773  widget_control, base2, sensitive = 1 
    6874 
     75  print, 'average clumping factor:', tclump 
     76 
    6977return 
    7078end 
Note: See TracChangeset for help on using the changeset viewer.