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 mksim(pbsbase,**d):
    print("making {}".format(pbsbase));
    myd = sd(lsp_d, **d);
    lsp=genlsp(**myd);
    #two color hack
    if test(myd, 'relative_phase'):
        relative_phases="phase {}".format(myd['relative_phase'])
    else:
        relative_phases=''
    lasers = '''
laser
wavelength 780e-7
spotsize   2.17e-4

laser
amplitude  0.36
wavelength 390e-7
{phase}
end
'''.format(phase=relative_phases);
    if test(myd, 'two_colors'):
        lsp = re.sub(r'type *19.*','type 85',lsp);
        lsp = re.sub(r' *; *\\lambda *spotsize','',lsp);
        lsp = re.sub('.*coefficients 7.80*e-03.*$', lasers,lsp,flags=re.MULTILINE);
    pbs=genpbs(pbsbase=pbsbase,domains=myd['domains']);
    pbs = re.sub("../scripts/autozipper","../../scripts/autozipper",pbs);
    pbs = re.sub("lsp-10-xy","lsp-10-xy-multilaser",pbs);
    output(lsp,pbs,pbsbase,
           dats=["sine700points.dat", myd['targetdat']],
           dir=pbsbase);

def genregions(**kw):
    if test(kw,'region_split'):
        regs = genregions_uniform([kw['region_split']],**kw);
    elif test(kw, 'region_splits'):
        regs = genregions_uniform(kw['region_splits'],**kw);
    else:
        regs = genregions_uniform([],**kw);
    return mkregion_str(regs);

def gendens(**kw):
    getkw = mk_getkw(kw, densdefaults);
    speciesl =  getkw('speciesl');
    outputfmt = "n_{}";
    if test(kw,'target_density'):
        Q = kw['target_density'];
        fracs = getkw('fracs');
        if test(kw, 'target_density_plainconst'):
            ret = {
                outputfmt.format(species) : plainconst_tmpl.format(data=Q*f)
                for species,f in zip(speciesl,fracs) };
            kw.update(ret);
            return kw;
        if type(Q) == tuple:
        #the reason for tuple and not a general iterable
        #is when we pass a single scale, for example, which
        #we scale by fracs
            pass;
        else:
            if hasattr(Q[0], '__call__'):
                Q = [  lambda x: frac*iQf(x)
                       for frac,iQf in zip(fracs,Q) ];
            else:
                Q = [ frac*iQ
                      for frac,iQ in zip(fracs,Q) ];
        if hasattr(Q[0],'__call__'):
            #invariant: txmin,txmax are in cm.
            x = np.linspace(kw['txmin'][0],kw['txmax'][1],20);
            Q[:] = [iQ(x) for iQ in Q];
        ret = {
            outputfmt.format(species) : densitypairs_tmpl.format(data = iQ)
            for species,iQ in zip(speciesl,Q) };
        kw.update(ret)
        return kw;
    else:
        kw['dens_dat'] = getkw('dens_dat');
        kw['dens_imul'] = getkw('dens_imul');
        kw['dens_type'] = getkw('dens_type');
        kw['fracs']     = getkw('fracs');
        def copy_kw(l):
            if type(kw[l]) != tuple:
                kw[l] = (kw[l],)*len(speciesl)
        copy_kw('dens_dat');
        copy_kw('dens_imul');
        copy_kw('dens_type');
        for i,species in enumerate(speciesl):
            kw['n_'+species] = densityfile_tmpl.format(
                targetdat = kw['dens_dat'][i],
                type = kw['dens_type'][i],
                imul = kw['dens_imul'][i],
                dmul = kw['fracs'][i]);
        return kw;
    pass;

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

def gengrids(**kw):
    getkw = mk_getkw(kw,grid_defs,prefer_passed = True);
    outgrids = ['','',''];
    simple_gridtmpl = '''
{dim}min             {min:e}
{dim}max             {max:e}
{dim}-cells          {cells}''';
    vargrid_tmpl = '''
{dim}min             {min:e}
{dim}max             {max:e}
{dim}-cells          {cells}
{dim}-intervals
 d{dim}-start {dxstart:e}
{intervals}
end'''
    vargrid_intv_tmpl = " length {length:e} for {N}";
    grids = []
    for dim in 'xyz':
        if getkw('{}cells'.format(dim)) == 0:
            grids.append("");
            continue;
        vgridlabel = '{}vargrid'.format(dim);
        if test(kw,vgridlabel):
            intv = '\n'.join([
                vargrid_intv_tmpl.format(length=l, N=N)
                for l,N in getkw(vgridlabel)]);
            if test(kw, 'd{}start'.format(dim)):
                ddimstart = getkw('d{}start'.format(dim));
            else:
                l = getkw(vgridlabel)[0][0];
                N = getkw(vgridlabel)[0][1];
                ddimstart = l/N;
            grid = vargrid_tmpl.format(
                dim=dim,
                min=getkw('{}min'.format(dim)),
                max=getkw('{}max'.format(dim)),
                cells=getkw('{}cells'.format(dim)),
                dxstart=ddimstart,
                intervals=intv)
        else:                    
            grid= simple_gridtmpl.format(
                dim=dim,
                min=getkw('{}min'.format(dim)),
                max=getkw('{}max'.format(dim)),
                cells=getkw('{}cells'.format(dim)));
        grids.append(grid);
    return grids;

def gentemp(**kw):
    getkw = mk_getkw(kw, tempdefaults);
    speciesl = getkw('speciesl');
    otherfuncs = '';
    funcnum = getkw('funcnum');
    if test(kw,'target_temps'):
        Q = kw['target_temps'];
        def process_temp(iq):
            if not iq: return None;
            if hasattr(iq, '__call__'):
                xres = getkw('dat_xres')
                x = np.linspace(kw['txmin'][0],kw['txmax'][1],xres);
                iq = iq(x);
                raise ValueError("Not implemented yet!");
            elif type(iq) is dict:
                _getkw = mk_getkw(iq, species_tempdefault);
                ret = densityfile_tmpl.format(
                    targetdat = _getkw('dat'),
                    type = _getkw('type'),
                    imul = _getkw('imul'),
                    dmul = _getkw('frac'));
            return ret;                
        ss = [ process_temp(iq) for iq in Q ];
    else:
        Q = [dict()]*len(speciesl);
        ss= [None]  *len(speciesl);
    for iq,species,s,e in zip(Q,speciesl,ss,getkw('thermal_energy')):
        cur = 'thermal_energy {}\n'.format(e);
        if s:
            otherfuncs += "function{}\n".format(funcnum);
            otherfuncs += s;
            # if 'energy_flags' in iq:
            #     energyflags = iq['energy_flags'];
            # else:
            #     energyflags = getkw('energy_flags');
            cur += 'spatial_function {}\n'.format(funcnum);
            # cur += 'energy_flags {}\n'.format(
            #     joinspace([
            #         1 if i else 0
            #         for i in getkw("dens_flags")]));
            funcnum += 1;
        kw['{}_thermalopts'.format(species)] = cur;
    if otherfuncs != '':
        if not test(kw, 'other_funcs'):
            kw['other_funcs'] = ''
        kw['other_funcs'] += otherfuncs;
    return kw;

 def format_region(region):
     if not test(region,"cells"):
         region['cells'] = ""
     else:
         region['cells'] = "cells = {}".format(region['cells'])
     if split_cells is not None and split_cells/region['domains'] < 4:
         print("warning: {} limits less than 4 cells thick".format(
             region['split'][0]));
     return region_tmpl.format(**region);

 def setspeciesv(fmt,defaultl,join=False,scale=False):
     l = fmt.format(species);
     if test(kw,l):
         v = kw[l];
         if scale: v=mt(v,getkw('ux'));
         if join:  v=joinspace(v);
     else:
         v = fmtd[defaultl];
     fmtd[l] = v;

def genlsp(**kw):
    def getkw(l,scale=None):
        if test(kw, l):
            if scale:
                return [scale*i for i in kw[l]];
            return kw[l];
        else:
            return defaults[l];
    E0 = np.sqrt(2*getkw('I')*1e4/(c*e0))*1e-5
    xmin,xmax, ymin,ymax = getkw('lim',scale=1e-4)
    fp = joinspace(getkw("fp"));
    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)*1e9,
        ((ymax-ymin)/ycells/c)*1e9)
    if timestep > couraunt:
        import sys
        sys.stderr.write("warning: timestep exceeds couraunt limit\n");
    targetdat = getkw('targetdat');
    dumpinterval=getkw('dumpinterval')*1e9;
    description=getkw('description');
    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'),
        domains=domains,totalt=totalt,
        timestep=timestep,
        targetdat=targetdat,
        dumpinterval=dumpinterval,
        description=description
    );
    return s;

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 genregions_uniform(subdivs,**kw):
    '''
    genregions in a uniform fashion

    subdivs is a list of tuples, each of which specify a range to divide
    along an axis. There are two options, the first is
       (d, split)
    in which d is a string representing the dimension, and split is
    the number of subdivions along that axis.
    The section option is
       (d, split, dmin, dmax) or (d, split, (dmin, dmax))
    where dmin is the minimum along the d axis and dmax is the max.

    The subdivisions are cartesian multiplied (or cartesian product).
    '''
    getkw = mk_getkw(kw, unireg_defaults);
    lims = mt(getkw('lim'),getkw('ux'));
    dims = 'xyz';
    out={}
    for subdiv in subdivs:
        if len(subdiv) == 3:
            subdiv = subdiv[:2] + subdiv[2];
        elif len(subdiv) == 2:
            i=dims.index(subdiv[0])
            subdiv = subdiv + tuple(lims[i*2:i*2+2])
        ax, split, mn, mx = subdiv;
        if ax not in out:
            out[ax]=[];
        out[ax].extend(list(
            np.linspace(mn,mx,split+1)));
    for i,dim in enumerate(dims):
        if dim not in out:
            out[dim] =lims[i*2:i*2+2];
    regs = [ dict(xmin=xmin,xmax=xmax,
                  ymin=ymin,ymax=ymax,
                  zmin=zmin,zmax=zmax,)
             for xmin,xmax in zip(out['x'],out['x'][1:])
             for ymin,ymax in zip(out['y'],out['y'][1:])
             for zmin,zmax in zip(out['z'],out['z'][1:]) ];
    for i,reg in enumerate(regs):
        reg['i'] = i+1;
    if test(kw,"domains_per_region"):
        for reg in regs:
            reg['domains'] = getkw("domains_per_region");
    else:
        ndom = getkw("domains")//len(regs);
        ex   = getkw("domains") %len(regs);
        for reg in regs:
            reg['domains'] = ndom;
        regs[-1]['domains'] += ex;
    for reg in regs:
        reg['split'] = "{}SPLIT".format(
            getkw("region_dom_split").upper());
    return regs;

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 highlight(ret, val,
              q=None, color='white', alpha=0.15, erase=False):
    '''
    Highlight a pc. Essentially a wrapper of plt.contour
    
    Arguments:
      ret   -- dict returned from pc.
      val   -- value to highlight
      q     -- quantity to highlight. If None, highlight ret's quantity
    
    Keyword Arguments:
      color -- color of highlight
      alpha -- alpha of highlight
      erase -- erases the highlights. Defaults to false (opposite of matplotlib!)
    
    Returns:
      ret but with stuff that plt.contour adds.
    '''
    ax = ret['axes'];
    x,y=ret['x'],ret['y'];
    if q is None:
        q = ret['q'];
    if test(ret,'flip') or test(ret,'rotate'):
        x,y=y,x;
    #elif q is not ret['q'] and test(ret,'flip'):
    if not test(ret, 'cbar'):
        ret['cbar'] = plt.colorbar(ret['pc']);
    cbar = ret['cbar'];
    if not test(ret, 'cts'):
        ret['cts'] = [];
    ct = ax.contour(x,y,q, [val],
                    colors=[color], alpha = alpha);
    ret['cts'].append(ct);
    if q is ret['q']:
        cbar.add_lines(ct,erase=erase);
    return ret;

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;

t = d['t'];
tbins = np.arange(ti,t[-1]+tstep,tstep);
#fucking c like loop shit mother fucker.
i=0;
Is=[];
for j,ct in enumerate(t):
    if ct > tbins[i]:
        Is.append(j);
        i+=1;
#do first
surf,ax,fig,bins,data = angular(d,**kw);
kw['fig'] = fig;
rmax=ax.get_rmax()
t=fig.text(tx,ty,'t = {:4.1f} fs'.format(tbins[0]*1e6),
           fontdict={'fontsize':22});
if not test(kw, 'phi'):
    kw['phi']='phi'

def animate(ii):
    j,i = ii;
    #plt.clf();
    #_,ax,_,_ = angular(d[:i],**kw);
    #ax.set_rmax(rmax);
    S,_,_ = np.histogram2d(
        data[0][:i],
        data[1][:i],
        bins=bins,
        weights=data[2][:i]);
    surf.set_array(S[::,:-1].ravel());
    #t=fig.text(tx,ty,'t = {:3.2f}e-4 ns'.format((tbins[j]-mt)*1e4),
    #       fontdict={'fontsize':22});

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,
#.........这里部分代码省略.........

def genoutlets(**kw):
    outlet_tmpl='''
;{label}
outlet
from {xmin:e}  {ymin:e} {zmin:e}
to   {xmax:e}  {ymax:e} {zmax:e}
phase_velocity 1.0
drive_model NONE''';
    laser10_tmpl='''
;laser
outlet
from {xmin:e}  {ymin:e} {zmin:e}
to   {xmax:e}  {ymax:e} {zmax:e}
phase_velocity 1.0
drive_model LASER
reference_point {fp}
components {components}
phases {phases}
temporal_function {lasertfunc}
analytic_function {laserafunc}
time_delay {time_delay}
'''
    retoutlets = ''
    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'):
        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);
        retoutlets += laser10_tmpl.format(
            fp = joinspace(scaletuple(lgetkw("fp"))),
            components = joinspace(lgetkw("components")),
            phases =  joinspace(lgetkw("phases")),
            lasertfunc = lgetkw('lasertfunc'),
            laserafunc = lgetkw('laserafunc'),
            time_delay = lgetkw('laser_time_delay'),
            **outlet_coords(outlet, l)
        );
    #outlet boundaries
    for side in [i for i in all_lims if i not in lset] :
        if laserkw[side[0]+'cells'] > 0:
            retoutlets += outlet_tmpl.format(
                label = side_label[side],
                **outlet_coords(side, laserkw));
    return retoutlets;