function prepare_gridgen_input(rotation, xoffset, yoffset, nrow, ncol, dx, dy, refels, varargin)
    % Refels is an array of structs (refinement elements). Each struct has a
    % 2d array of points [x0,y0;x1,y1;...], and a property ref_level that
    % gives the level of refinement. If the array of points has a single
    % entry, it will be assumed to be a point. If the last point has the same
    % coordinates as the first point, it is assumed to be a polygon.
    % Otherwise, it is assumed to be a polyline.

    mk_batch_file=0;
    for ii=1:length(varargin)
        if ischar(varargin{ii})
            switch lower(varargin{ii})
                case{'make_batch'}
                    mk_batch_file=1;
            end
        end
    end

    fn_out = 'qtg.dfn'; % Output filename(s)

    [~, bn_out, ~] = fileparts(fn_out);

    fn_in_1 = 'qtg1.dfn'; % The main input file to be generated
    fn_in_2 = 'qtg2.dfn'; % The generated input file to convert grid to shp file

    fn_bat = 'run_gridgen.bat';

    % TODO: Assert no spaces in filenames

    main_fmt = ['begin modflow_grid basegrid\n' ...
          '  rotation_angle = %f\n' ...
          '  x_offset = %d\n' ...
          '  y_offset = %d\n' ...
          '  length_unit = meter\n' ...
          '  nlay = 1\n' ...
          '  nrow = %d\n' ...
          '  ncol = %d\n' ...
          '  delr = constant %f\n' ...
          '  delc = constant %f\n' ...
          '  top = constant 0.000000\n' ...
          '  bottom layer 1 = constant -100.000000\n' ...
          'end modflow_grid\n' ...
          '\n' ...
          'begin quadtree_builder qtgbuild\n' ...
          '  modflow_grid = basegrid\n' ...
          '%s\n' ...
          '  smoothing = full\n' ...
          '  smoothing_level_horizontal = 1\n' ...
          '  smoothing_level_vertical = 1\n' ...
          '  top layer 1 = replicate basegrid\n' ...
          '  bottom layer 1 = replicate basegrid\n' ...
          '  grid_definition_file = %s\n' ...
          'end quadtree_builder\n'];

    ref_layers_str = '';

    if ~isempty(refels)
        ref_layers_str = [ref_layers_str '  refinement_features layer 1 ='];
    end
    for i=1:length(refels)
        if size(refels(i).xy, 1) == 1
            % This is a point
            refels(i).type = 1;
            refels(i).typename = 'point';
        elseif all(refels(i).xy(1,:) == refels(i).xy(end,:))
            % This is a polygon
            refels(i).type = 5;
            refels(i).typename = 'polygon';
        else
            % Assumed polyline
            refels(i).type = 3;
            refels(i).typename = 'polyline';
        end
        tmp_str = sprintf(' %s%d', refels(i).typename, i);
        ref_layers_str = [ref_layers_str tmp_str];
    end

    main_str = sprintf(main_fmt, rotation, xoffset, yoffset, ...
                       nrow, ncol, dx, dy, ref_layers_str, fn_out);

    % Generate and append point/polyline/polygon shapefiles
    refel_fmt = ['\nbegin refinement_features %s%d\n' ...
                    '  shapefile = %s%d\n' ...
                    '  feature_type = %s\n' ...
                    '  refinement_level = %d\n' ...
                    'end refinement_features\n'];

    for i=1:length(refels)
        el = refels(i);
        xy = {el.xy};
        gridgen_type = el.typename;

        if el.type == 1
            % Points go directly, not in a cell array
            xy = el.xy;
        end
        if el.type == 3
            % Gridgen type is call "line", not polyline
            gridgen_type = 'line';
        end

        cur_str = sprintf(refel_fmt, el.typename, i, ...
                                     el.typename, i, ...
                                     gridgen_type, el.ref_level);

        main_str = [main_str cur_str];
        shapewrite([el.typename num2str(i) '.shp'], el.typename, xy);
    end

    % And write everything to file
    fid = fopen(fn_in_1, 'wt');
    fprintf(fid, main_str);
    fclose(fid);

    to_shp_fmt = ['load %s\n' ...
        '\n' ...
        'begin grid_to_shapefile qtg-to-shapefile\n' ...
        '  grid = %s\n' ...
        '  shapefile = %s\n' ...
        '  feature_type = polygon\n' ...
        'end grid_to_shapefile\n' ...
        '\n' ...
        'begin grid_to_usgdata qtg-to-crs\n' ...
        '  grid = %s\n' ...
        '  usg_data_prefix = %s\n' ...
        'end grid_to_usgdata\n'];

    to_shp_str = sprintf(to_shp_fmt, fn_out, bn_out, bn_out, bn_out, bn_out);

    fid = fopen(fn_in_2, 'wt');
    fprintf(fid, to_shp_str);
    fclose(fid);

    if mk_batch_file
        bat_fmt = [':: construct grid\n' ...
            'gridgen_x64.exe qtgbuild %s\n' ...
            '\n' ...
            ':: write grid shapefile\n' ...
            'gridgen_x64.exe qtg-to-shapefile %s\n' ...
            '\n' ...
            ':: write grid nod-file\n' ...
            'gridgen_x64.exe qtg-to-crs %s\n'];
        
        fid = fopen(fn_bat, 'wt');
        fprintf(fid, bat_fmt, fn_in_1, fn_in_2, fn_in_2);
        fclose(fid);
    end

end

