def gendats(di,
            w0=w0*1e2,
            width=0.46e-4,
            L=0.043e-4,
            N0=1.08e22,
            depth=0.086e-4,
            mindensity=1e18,
            dat_xres=None,
            new=False):
    if new:
        mkpinprick = mk45_pinprick_plasma
    else:
        mkpinprick = mk45_pinprick_plasma_old
    targ_plasma, targ_neutral = mkpinprick(
        dim = [i*1e-4 for i in d['tlim']],
        N0  = N0,
        laser_radius = w0,
        width = width,
        L = L,# 43nm
        depth = depth, #chosen arbitrarily
        mindensity=mindensity);
    if not dat_xres:
        dat_xres = di['res'][0]+1;
    print("making targets for {}".format(di['pbsbase']));
    dd = sd(di, f_2D = targ_plasma, dat_xres = dat_xres);
    dat = gendat(**dd);
    savetxt(
        "{}/{}".format(di['pbsbase'],'target_plasma.dat'),
        dat);
    dd = sd(d, f_2D = targ_neutral, dat_xres = dat_xres);
    dat = gendat(**dd);
    savetxt(
        "{}/{}".format(di['pbsbase'],'target_neutral.dat'),
        dat);

def manual_genboundary(bspec,**kw):
    kwp = sd(kw,**bspec);
    getkw = mk_getkw(kwp,outletdefaults,prefer_passed = True);
    btype = getkw("type");
    lims  = getkw('lim');
    di = dict();
    for dim,lim in zip(all_lims,lims):
        di[dim] = lim;
    if btype == 'outlet':
        model = getkw('model');
        di['label']= getkw('label');
        di['phase_velocity'] = getkw('phase_velocity');
        if model == 'none':
            ret = outlet_none_tmpl.format(**di);
        elif model == 'potential':
            di = sd(manbounds_defs,**di);
            if not test(di,'connection_rank'):
                di['connection_rank'] ='';
            else:
                di['connection_rank'] = 'connection_rank {}'.format(
                    getkw('connection_rank'));
            if not test(kwp,'voltage_measurement'):
                di['voltage_measurement'] = '';
            else:
                raise NotImplementedError("haven't got to this yet...");
            di['circuit'] = di['circuit'];
            ret = outlet_pot_tmpl.format(**di);
        elif model == 'laser':
            raise NotImplementedError("need to implement the laser...");
        else:
            raise ValueError("unknown outlet model {}".format(model));
    else:
        raise NotImplementedError("yeah...");
    return ret;

def genregions(**kw):
    def getkw(l, scale=None):
        if test(kw, l):
            ret = kw[l]
        else:
            ret = region_defaults[l]
        if scale:
            return [scale * i for i in ret]
        return ret

    regsplit_dir, subdivs = getkw("region_split")

    nonsplits = [x for x in ["x", "y", "z"] if x != regsplit_dir]
    limkw = [x + lims for x in ["x", "y", "z"] for lims in ["min", "max"]]
    lims = {k: lim for k, lim in zip(limkw, getkw("lims"))}

    total_doms = getkw("domains")
    doms = [total_doms // subdivs for i in range(subdivs)]
    doms[-1] += total_doms % subdivs
    lmn, lmx = regsplit_dir + "min", regsplit_dir + "max"
    mn, mx = lims[lmn], lims[lmx]
    edges = [mn + i * (mx - mn) / subdivs for i in range(subdivs)] + [mx]
    mins = edges[:-1]
    maxs = edges[1:]
    reg = sd(lims, split=getkw("split_dir").upper + "SPLIT")
    regions = [
        sd(reg, **{lmn: mn, lmx: mx, "i": i, "domains": di}) for i, (mn, mx, di) in enumerate(zip(mins, maxs, doms))
    ]
    return mkregion_str(regions)

def gendat3d(
        di,
        w0=w0*1e2,
        width=0.46e-4,
        L=0.043e-4,
        N0=1.08e22,
        mindensity=1e18,
        dat_xres=None,
        dat_zres=4,#works for no z variation
        fmt='%.4e'
):
    targ_neutral = mk45_pinprick_neutral3d(
        dim = [i*1e-4 for i in d['tlim']],
        N0  = N0,
        laser_radius = w0,
        width = width,
        L = L,# 43nm
        mindensity=mindensity);
    if not dat_xres:
        dat_xres = di['res'][0]+1;
    print("making targets for {}".format(di['pbsbase']));
    dd = sd(di, f_3D = targ_neutral, dat_xres = dat_xres);
    dat = gendat(dat_zres=dat_zres,datfmt=fmt,**dd);
    savetxt(
        "{}/{}".format(di['pbsbase'],di['dens_dat']),
        dat);

def genonescale(**kw):
    getkw = mk_getkw(kw, onescale_defaults)
    slen = getkw("solid_len")
    xlen = getkw("xlen")
    kw1 = sd(kw, tlim=(0.0, xlen) + (0.0, 0.0, 0.0, 0.0), sdim=(xlen - slen, xlen) + (0.0, 0.0, 0.0, 0.0))
    kw1["f_1D"] = genf(**kw1)
    return gentargetdat(**kw1)

def mksim(pbsbase, **d):
    print("making {}".format(pbsbase))
    myd = sd(lsp_d, **d)
    lsp = genlsp(**myd)
    pbs = genpbs(pbsbase=pbsbase)
    pbs = re.sub("../scripts/autozipper", "../../scripts/autozipper", pbs)
    output(lsp, pbs, pbsbase, dats=["sine700points.dat", myd["targetdat"]], dir=pbsbase)

def firsthash_new(frame,**kw):
    kw['new']=True;
    kw['dupes']=None;
    hashes = genhash(frame,**kw);
    uni,counts = np.unique(hashes,return_counts=True);
    d=sd(kw,dupes=uni[counts>1],removedupes=True);
    dupei = np.in1d(hashes, d['dupes'])
    hashes[dupei] = -1
    return hashes, retd;

def mkscale_sim(**d):
    datf = "water-{}.dat".format(d["dat"])
    d = sd(scale_sim, **d)
    d["targetdat"] = datf
    d["description"] = "scale length sim with {}".format(d["dat"])
    if not test(d, "name"):
        name = "longl_l={}_{}".format(d["dat"], d["I"])
    else:
        name = d["name"]
    mksim(name, **d)

def genregions(**kw):
    def getkw(l,scale=None):
        if test(kw, l):
            ret = kw[l]
        else:
            ret = region_defaults[l];
        if scale:
            return [scale*i for i in ret];
        return ret;
    regsplit_dir, subdivs = getkw('region_split');
    
    nonsplits = [x for x in ['x','y','z']
                 if x != regsplit_dir ];
    limkw = [x+lims
             for x in ['x','y','z']
             for lims in ['min','max']];
    lims = {k:lim for k,lim in zip(limkw,getkw('lim',scale=1e-4))};
    
    total_doms = getkw('domains');
    doms =  [total_doms//subdivs for i in range(subdivs)];
    doms[-1] += total_doms % subdivs;
    lmn,lmx = regsplit_dir+'min', regsplit_dir+'max';
    mn,mx = lims[lmn],lims[lmx]
    edges = [mn+i*(mx-mn)/subdivs
             for i in range(subdivs)] + [mx];
    split_cells = getkw(getkw("region_dom_split")+"cells");
    
    mins = edges[:-1];
    maxs = edges[1:];
    if test(kw,"xcells") and test(kw,"ycells") and test(kw,"zcells"):
        zcells_per_region = [kw['zcells']//subdivs 
                             for i in range(subdivs)]
        zcells_per_region[-1] += kw['zcells'] % subdivs;
        
        cellses = [ kw['xcells']*kw['ycells']*zc
                    for i,zc in enumerate(zcells_per_region)];
    else:
        cellses = [None for i in range(subdivs)];
    reg = sd(lims,split=getkw("region_dom_split").upper()+"SPLIT")
    regions = [ sd(reg,**{lmn:mn,lmx:mx,'i':i+1,'domains':di,'cells':cells})
                for i,(mn,mx,di,cells) in enumerate(zip(mins,maxs,doms,cellses)) ];
    return mkregion_str(regions, split_cells=split_cells);

 def mktarg(di):
     dd = sd(
         di,
         f_2D = mk45(
             dim = [i*1e-4 for i in di['tlim']],
             N0    = 1.0804e22,
             width = 0.46e-4,
             dropcorners='round'));
     dat = gendat(**dd);
     savetxt(
         "{}/{}".format(di['pbsbase'],di['dens_dat']),
         dat);

def genconductor_boundaries(**kw):
    getkw=mk_getkw(kw, condb_defaults);
    conductorss='';
    for I,conductor in enumerate(getkw('conductors')):
        cd = sd(condb_objdef, **conductor);
        if test(cd,'from') and test(cd,'to'):
            cd['xmin'],cd['ymin'],cd['zmin'] = cd['from']
            cd['xmax'],cd['ymax'],cd['zmax'] = cd['to']
            pass;
        else:
            outlet = cd['outlet'];
            if outlet not in all_lims:
                raise ValueError('Unknown outlet "{}"'.format(outlet));
            coords = outlet_coords(outlet,kw);
            cd['width']*=1e-4;
            cd['start']*=1e-4;
            sign = lambda outlet: 1.0 if outlet[-2:] == 'ax' else -1.0
            coords[outlet] += sign(outlet)*(cd['width'] + cd['start']);
            coords[otherside(outlet)] += sign(outlet)*cd['start'];
        conductorss += condb_tmpl.format(
            i=I+1,
            **sd(cd,**coords));
    return conductorss;

def gendatrot(
        di,
        w0=w0*1e2,
        width=0.46e-4,
        L=0.043e-4,
        N0=1.08e22,
        mindensity=1e18,
        spotz_width=20e-4,
        yres = 75,
        fmt='%.4e',):
    print("warning: this will only work for rot3d");
    targ_neutral = mk0_pinprick_neutral3d(
        N0  = N0,
        spotz_width=spotz_width,
        laser_radius = w0,
        width = width,
        L = L,# 43nm
        mindensity=mindensity);
    #manual hacks to save space
    #generate cell sized samples around the corners.
    # x's edge
    def getsamples(point,lims,res,ndxs = 3):
        out=np.linspace(lims[0]*1e-4,lims[1]*1e-4,res+1);
        dx = out[1] - out[0];
        out = out[ out >= point - ndxs*dx ];
        out = out[ out <= point + ndxs*dx ];
        return list(out);
    def axissamples(half_width, tlim, lim, res,ndxs=3):
        p = [ tlim[0] ] + getsamples(-half_width,lim,res,ndxs=ndxs)
        p+= [0.0] + getsamples(half_width,lim,res,ndxs=ndxs)
        p+= [ tlim[1] ];
        return np.array(p);
    x = axissamples(
        w0, di['tlim'][0:2], di['lim'][0:2], di['res'][0]);
    y = np.linspace(
        di['tlim'][2]*1e-4,di['tlim'][3]*1e-4,yres+1);
    z = axissamples(
        spotz_width/2.0,
        di['tlim'][4:6], di['lim'][4:6], di['res'][2]);
    print("making targets for {}".format(di['pbsbase']));
    X,Y,Z=np.meshgrid(x,y,z,indexing='ij');
    Q=targ_neutral(X,Y,Z)
    dd = sd(di, data3D =(x,y,z,Q), datfmt=fmt);
    dat = gendat(**dd);
    savetxt(
        "{}/{}".format(di['pbsbase'],di['dens_dat']),
        dat);

def genconductors(**kw):
    getkw=mk_getkw(kw, condb_defaults);
    conductorss='';
    for I,conductor in enumerate(getkw('conductors')):
        cd = sd(condb_objdef, **conductor);
        coords=dict();
        if test(cd,'outlet'):
            outlet = cd['outlet'];
            if outlet not in all_lims:
                raise ValueError('Unknown outlet "{}"'.format(outlet));
            coords = outlet_coords(outlet,kw);
            cd['width']*=1e-4;
            cd['start']*=1e-4;
            if outlet[-2:] == 'ax':
                sign = 1.0;
            else:
                sign =-1.0
            cd['type']= 'BLOCK';
            coords[outlet] += sign*(cd['width'] + cd['start']);
            coords[otherside(outlet)] += sign*cd['start'];
            cd['from'] = (coords['xmin'],coords['ymin'],coords['zmin']);
            cd['to']   = (coords['xmax'],coords['ymax'],coords['zmax']);
        conductorss += condf_tmpl.format(
            i=I+1,
            xf=cd['from'][0],yf=cd['from'][1],zf=cd['from'][2],
            **cd);
        def mk_to(cd):
            '''generate the to's'''
            return ''.join([
                condt_tmpl.format(xt=xt,yt=yt,zt=zt)
                for xt,yt,zt in cd['to'] ]);
        
        if cd['type'] == 'BLOCK':
            if type(cd['to']) == list: cd['to'] = cd['to'][0]
            conductorss += condt_tmpl.format(
                xt=cd['to'][0],yt=cd['to'][1],zt=cd['to'][2]);
        elif cd['type'] == 'PARALLELPIPED':
            conductorss += mk_to(cd);
        elif cd['type'] == 'TRILATERAL':
            conductorss += mk_to(cd);
            conductorss += 'sweep_direction {sweep_direction}\n'.format(
                **cd);
        else:
            raise ValueError(
                "Unknown conductor type '{}'".format(cd['type']));
    return conductorss;

def gendatclean(
        di,
        width=0.46e-4,
        N0=1.08e22,
        fmt='%.4e',):
    targ_neutral = mk45_clean_neutral2d(
        width=width,
        N0  = N0,);
    print("making targets for {}".format(di['pbsbase']));
    dd = sd(di,
            f_2D = targ_neutral,
            dat_xres=di['res'][0]+1,
            dat_yres=di['res'][0]+1);
    dat = gendat(datfmt=fmt,**dd);
    savetxt(
        "{}/{}".format(di['pbsbase'],di['dens_dat']),
        dat);

        contour_quantities=('RhoN10', 'RhoN10'),
    ),
    #pmovies
    no_pmovies=True,
    #particle dumps
    dump_particle=True,
    particle_dump_interval_ns=1e-15,
);
d2 = sd(d,
        lim = (-25,25,
               -25,25,
               0,0),
        tlim = (-20,20,
                -20,20,
                0,0),
        res =(1600,
              1600,
              0),
        pbsbase='prexp2',
        region_split=('y',3),
        domains=96,
        totaltime=400e-15,
        particle_dump_interval_ns=1e-15,);
gensim(**d);
gensim(**d2);

def mktarg(di, No = 3.34e22, ro = 5e-4, L = 1e-4):
    tw = di['tlim'][1] - di['tlim'][0]
    dat_xres= int(tw/(di['lim'][1]-di['lim'][0])*di['res'][0]);
    dd = sd(
        di,

def genlsp(**kw):
    def getkw(l, scale=None):
        if test(kw, l):
            ret = kw[l]
        else:
            ret = defaults[l]
        if scale:
            return [scale * i for i in ret]
        return ret

    intensity = getkw("I")
    E0 = np.sqrt(2 * getkw("I") * 1e4 / (c * e0)) * 1e-5
    xmin, xmax, ymin, ymax = getkw("lim", scale=1e-4)
    fp = joinspace(getkw("fp", scale=1e-4))
    components = joinspace(getkw("components"))
    phases = joinspace(getkw("phases"))
    l = getkw("l") * 100.0
    if test(kw, "resd"):
        xres, yres = getkw("resd")
        xcells = (xmax - xmin) / (l / xres)
        ycells = (ymax - ymin) / (l / yres)
    else:
        xcells, ycells = getkw("res")
    w0 = getkw("w") * 100.0
    T = getkw("T") * 1e9
    targ_xmin, targ_xmax, targ_ymin, targ_ymax = getkw("tlim", scale=1e-4)
    domains = getkw("domains")
    # we have that na~l/(pi*w), and the f-number~1/2na, thus
    # f-number ~ pi*w/2l
    fnum = np.pi * w0 / 2 / l
    totalt = getkw("totaltime") * 1e9
    timestep = getkw("timestep") * 1e9
    couraunt = min(((xmax - xmin) / xcells / c_cgs) * 1e9, ((ymax - ymin) / ycells / c_cgs) * 1e9)
    if timestep > couraunt:
        import sys

        sys.stderr.write("warning: timestep exceeds couraunt limit\n")
    pexts = genpext(**sd(kw, species=getkw("pext_species")))
    targetdat = getkw("targetdat")
    dumpinterval = getkw("dumpinterval") * 1e9
    description = getkw("description")
    restart = getkw("restart")
    if not test(kw, "no_pmovies"):
        pmovies = """
particle_movie_interval_ns {dumpinterval}
particle_movie_components Q X Y VX VY XI YI
""".format(
            dumpinterval=dumpinterval
        )
    else:
        pmovies = ""
    restarts = "maximum_restart_dump_time {}".format(restart) if restart else ""
    with open("hotwater2d_tmpl.lsp") as f:
        s = f.read()
    s = s.format(
        xmin=xmin,
        xmax=xmax,
        ymin=ymin,
        ymax=ymax,
        xcells=xcells,
        ycells=ycells,
        l=l,
        w0=w0,
        E0=E0,
        fnum=fnum,
        targ_xmin=targ_xmin,
        targ_xmax=targ_xmax,
        targ_ymin=targ_ymin,
        targ_ymax=targ_ymax,
        fp=fp,
        pulse=T,
        components=components,
        phases=phases,
        intensity=getkw("I"),
        pmovies=pmovies,
        pexts=pexts,
        domains=domains,
        totalt=totalt,
        timestep=timestep,
        targetdat=targetdat,
        dumpinterval=dumpinterval,
        description=description,
        restarts=restarts,
    )
    return s

    "res": (1400, 1600),
    "tlim": (-27.5, 0, -15, 15),
    "fp": (0, 0, 0),
    "domains": 48,
    "totaltime": 300e-15,
    "timestep": 4e-17,
    "components": (0, 1, 0),
    "phases": (0, 0, 0),
    "targetdat": "watercolumn.dat",
    "dumpinterval": 2e-16,
    "description": "Hotwater in 2d",
    "pext_species": (10,),
    "restart": None,
}

pext_defaults = sd(defaults, species=(10,), start_time=0, stop_time=1)


def genpext(**kw):
    def getkw(l, scale=None):
        if test(kw, l):
            ret = kw[l]
        else:
            ret = pext_defaults[l]
        if scale:
            return [scale * i for i in ret]
        return ret

    tmpl = """
;
extract{i}

domains=1000;
region_split=('z',68);
pbsbase="flashic_3d";
defpbs = dict(
    pbsbase=pbsbase,
    pbsname=pbsbase+"_oakley",
    domains=domains,
    cluster='oakley',
    autozipper=False,
    queue=None,
    ppn=None,);
pbses=[
    defpbs,
    sd(
        defpbs,
        pbsname=pbsbase+"_garnet_debug",
        cluster='garnet',
        queue='debug'),
    sd(
        defpbs,
        pbsname=pbsbase+"_garnet",
        cluster='garnet',
        queue='standard_lw'),
    sd(
        defpbs,
        pbsname=pbsbase+"_garnet_short",
        cluster='garnet',
        queue='standard_lw',
        walltime=48),];
gensim(
    l = 0.8e-6,

def genobjects(**kw):
    getkw=mk_getkw(kw, condb_defaults);
    ux = getkw('ux');
    objss='';
    for I,objspec in enumerate(getkw('objects')):
        coords=dict();
        #objd = sd(condb_defaults, **kw);
        objd = sd(condb_objdef, **kw);
        objd = sd(objd, **objspec);
        objd['i'] = I + 1;
        if test(objspec,'outlet'):
            outlet = objd['outlet'];
            if outlet not in all_lims:
                raise ValueError('Unknown outlet "{}"'.format(outlet));
            coords = outlet_coords(outlet,kw);
            objd['width']*=ux;
            objd['start']*=ux;
            if outlet[-2:] == 'ax':
                sign = 1.0;
            else:
                sign =-1.0
            objd['type']= 'BLOCK';
            coords[outlet] += sign*(objd['width'] + objd['start']);
            coords[otherside(outlet)] += sign*objd['start'];
            objd['crossstart']*=ux;
            otherdims = [i for i in 'xyz' if i != outlet[:-3]];
            for dim in otherdims:
                if getkw('{}cells'.format(dim)) > 0:
                    coords['{}min'.format(dim)] += objd['crossstart'];
                    coords['{}max'.format(dim)] -= objd['crossstart'];
            objd['from'] = (coords['xmin'],coords['ymin'],coords['zmin']);
            objd['to']   = (coords['xmax'],coords['ymax'],coords['zmax']);
            #objss += condf_tmpl.format(
            #    i=I+1,
            #    condon=objd['condon'],
            #    xf=objd['from'][0],yf=objd['from'][1],zf=objd['from'][2],
            #    **objd);
        def mk_to(objd):
            '''generate the to's'''
            return ''.join([
                condt_tmpl.format(xt=xt,yt=yt,zt=zt)
                for xt,yt,zt in objd['to'] ]);
        frms = ['BLOCK', 'PARALLELPIPED', 'TRILATERAL'];
        if objd['type'] == 'SOLID':
            objss += obj_solid_tmpl.format(**objd);
        elif objd['type'] in frms:
            objss += condf_tmpl.format(
                xf=objd['from'][0],yf=objd['from'][1],zf=objd['from'][2],
                **objd);
            objd = sd(condb_objdef, **objd);
            if objd['type'] == 'BLOCK':
                if type(objd['to']) != list:
                    objd['to'] = [objd['to']];
            objss += mk_to(objd);
            if objd['type'] == 'TRILATERAL':
                objss += 'sweep_direction {sweep_direction}\n'.format(
                    **objd);
        elif objd['type'] == 'CONE':
            objd = sd(obj_cone_defs, **objd);
            coordtokeysl(objd, 'base', '{}b');
            coordtokeysl(objd, 'apex', '{}ap');
            coordtokeysl(objd, 'edge', '{}ed');
            objss += obj_cone_tmpl.format(**objd);
        elif objd['type'] == 'CYLINDER' or objd['type'] == 'TORUS':
            objd = sd(obj_cyl_defs, **objd);
            objd['axis'],objd['pitch'] = objd['axis_pitch'];
            objd['azaxis'],objd['azpitch'] = objd['azimuthal_axis_pitch'];
            objd['start_angle'],objd['sweep_angle'] = objd['azimuth_range'];
            if objd['type'] == 'CYLINDER':
                coordtokeysl(objd, 'base', '{}b');
                objss += obj_cyl_tmpl.format(**objd);
            elif objd['type'] == 'TORUS':
                objd = sd(obj_torus_defs,**objd);
                coordtokeysl(objd, 'center', '{}c');
                objss += obj_torus_tmpl.format(**objd);
            else:
                raise Exception("WTFdiosjf03y2q8qencdq");
        else:
            raise ValueError(
                "Unknown object type '{}'".format(objd['type']));
    pass #end for
    return objss;

def genboundaries(**kw):
    retboundaries = ''
    if test(kw,'manual_boundaries'):
        return '\n'.join([
            manual_genboundary(bspec,**kw)
            for bspec in kw['manual_boundaries'] ]);
    laserkw =  kw['laseroutlet'] if test(kw, 'laseroutlet') else dict();
    laserkw = sd(kw, **laserkw);
    getkw = mk_getkw(laserkw,outletdefaults,prefer_passed = True);
    lset = set();
    side_label = dict(
        xmin='front',
        xmax='back',
        ymin='left',
        ymax='right',
        zmin='bottom',
        zmax='top');
    ls = [];
    if not test(laserkw, 'multilaser') and not (test(laserkw,'nolaser') or test(laserkw,'nolaser_outlet')):
        ls = [ sd(laserkw, outlet='xmin') ];
    elif test(laserkw,'multilaser'):
        ls = kw['multilaser'];
        print("experimental multi-lasers");
        print("if you put more than one laser on a single outlet,");
        print("be sure to be using my modifications to lsp for it.");
        ls = kw['multilaser'];
    #lasers
    for l in ls:
        l = sd(laserkw, **l);
        lgetkw = mk_getkw(l, laserdefaults, prefer_passed = True);
        outlet = lgetkw('outlet');
        if outlet not in all_lims:
            raise ValueError('Unknown outlet "{}"'.format(outlet));
        lset.add(outlet);
        retboundaries += laser10_tmpl.format(
            fp = joinspace(mt(lgetkw("fp"),getkw('ux'))),
            components = joinspace(lgetkw("components")),
            phase_velocity = lgetkw('phase_velocity'),
            phases =  joinspace(lgetkw("phases")),
            lasertfunc = lgetkw('lasertfunc'),
            laserafunc = lgetkw('laserafunc'),
            time_delay = lgetkw('laser_time_delay'),
            **outlet_coords(outlet, l)
        );
    just_outlets = [i for i in all_lims if i not in lset];
    if test(kw,'freespace'):
        getkwfr = mk_getkw(
            sd(kw,**kw['freespace']),frsp_defs,prefer_passed=True);
        di=dict();
        di['model_type'] = getkwfr('model_type');
        if di['model_type'] not in ["WAVEABC","UNIAXIAL","CFSPML"]:
            raise ValueError("unrecognized model_type for freespace");
        if di['model_type']=="WAVEABC":
            di['num_of_cells'] = "";
        else:
            di['num_of_cells'] = "number_of_cells {}".format(
                getkwfr('num_of_cells'));
        keeps = getkwfr('keep_outlets');
        if keeps is None: keeps = [];
        just_outlets = [ i for i in just_outlets
                         if i in keeps ];
        if test(kw['freespace'],'frlim'):
            for dim,lim in zip(getkwfr('frlim'),all_lims):
                di[lim] = dim;
        else:
            dx = getkwfr('freesp_delta');
            keepset = lset.union(just_outlets);
            for lim in all_lims:
                di[lim] = getkw(lim);
                if lim in keepset:
                    if 'min' in lim:
                        di[lim] -= dx;
                    else:
                        di[lim] += dx;
        di['refp'] = joinspace(getkwfr('freesp_refp'));
        di['label']= getkwfr('freesp_label');
        retboundaries += freespace_tmpl.format(**di);
    #outlet boundaries which and are not removed by freespace
    for side in just_outlets:
        if laserkw[side[0]+'cells'] > 0:
            retboundaries += outlet_none_tmpl.format(
                label = side_label[side],
                phase_velocity=getkw('phase_velocity'),
                **outlet_coords(side, laserkw));
    pwbtmpl='''
planewave
from {xmin:e}  {ymin:e} {zmin:e}
to   {xmax:e}  {ymax:e} {zmax:e}
reference_point {refp}
polar_angle {polar}          
azimuthal_angle  {azimuth}
rotation_angle {rotation}
frequency {freq}
temporal_function {pwfunc}
'''
    pwbdefs = dict(
        polar=90,
        azimuth=45,
        rotation=-45,
        freq=1e3,
#.........这里部分代码省略.........