/* Steven Andrews, 1/16/03 */
/* Library of run-time interpretor commands for Smoldyn program. */
/* See documentation called Smoldyn doc */
/* Copyright 2003-2005 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "Rn.h"
#include "RnSort.h"
#include "Zn.h"
#include "random.h"
#include "string2.h"
#include "smollib2.h"
#include "opengl2.h"
#include "SimCommand.h"
#include "Geometry.h"

#include "opengl2.h"

int cmdstop(simptr sim,cmdptr cmd,char *line2);
int cmdpause(simptr sim,cmdptr cmd,char *line2);
int cmdbeep(simptr sim,cmdptr cmd,char *line2);
int cmdkeypress(simptr sim,cmdptr cmd,char *line2);

int cmdoverwrite(simptr sim,cmdptr cmd,char *line2);
int cmdincrementfile(simptr sim,cmdptr cmd,char *line2);

int cmdifno(simptr sim,cmdptr cmd,char *line2);
int cmdifless(simptr sim,cmdptr cmd,char *line2);
int cmdifmore(simptr sim,cmdptr cmd,char *line2);

int cmdpointsource(simptr sim,cmdptr cmd,char *line2);
int cmdvolumesource(simptr sim,cmdptr cmd,char *line2);
int cmdkillmol(simptr sim,cmdptr cmd,char *line2);
int cmdkillmolprob(simptr sim,cmdptr cmd,char *line2);
int cmdkillmolinsphere(simptr sim,cmdptr cmd,char *line2);
int cmdequilmol(simptr sim,cmdptr cmd,char *line2);
int cmdreplacexyzmol(simptr sim,cmdptr cmd,char *line2);
int cmdreplacevolmol(simptr sim,cmdptr cmd,char *line2);
int cmdmodulatemol(simptr sim,cmdptr cmd,char *line2);
int cmdreact1(simptr sim,cmdptr cmd,char *line2);
int cmdsetrateint(simptr sim,cmdptr cmd,char *line2);
int cmdsettimestep(simptr sim,cmdptr cmd,char *line2);
int cmdexcludebox(simptr sim,cmdptr cmd,char *line2);
int cmdexcludesphere(simptr sim,cmdptr cmd,char *line2);
int cmdincludeecoli(simptr sim,cmdptr cmd,char *line2);

int cmdecho(simptr sim,cmdptr cmd,char *line2);
int cmdmolcount(simptr sim,cmdptr cmd,char *line2);
int cmdlistmols(simptr sim,cmdptr cmd,char *line2);
int cmdlistmols2(simptr sim,cmdptr cmd,char *line2);
int cmdlistmols3(simptr sim,cmdptr cmd,char *line2);
int cmdmolpos(simptr sim,cmdptr cmd,char *line2);
int cmdmolmoments(simptr sim,cmdptr cmd,char *line2);
int cmdsavesim(simptr sim,cmdptr cmd,char *line2);
int cmdmeansqrdisp(simptr sim,cmdptr cmd,char *line2);

int insideecoli(double *pos,double *ofst,double rad,double length);
void putinecoli(double *pos,double *ofst,double rad,double length);
int molinpanels(simptr sim,int ll,int m,int s,char pshape);


/* ************ command processor ************ */

int docommand(void *cmdfnarg,cmdptr cmd,char *line) {
	simptr sim;
	char word[STRCHAR],*line2;
	int itct;

	if(!cmdfnarg) return 0;
	sim=(simptr) cmdfnarg;
	if(!line) return 0;
	itct=sscanf(line,"%s",word);
	if(!itct) return 0;
	line2=strnword(line,2);
	if(!strcmp(word,"stop")) return cmdstop(sim,cmd,line2);
	else if(!strcmp(word,"pause")) return cmdpause(sim,cmd,line2);
	else if(!strcmp(word,"beep")) return cmdbeep(sim,cmd,line2);
	else if(!strcmp(word,"keypress")) return cmdkeypress(sim,cmd,line2);

	else if(!strcmp(word,"overwrite")) return cmdoverwrite(sim,cmd,line2);
	else if(!strcmp(word,"incrementfile")) return cmdincrementfile(sim,cmd,line2);

	else if(!strcmp(word,"ifno")) return cmdifno(sim,cmd,line2);
	else if(!strcmp(word,"ifless")) return cmdifless(sim,cmd,line2);
	else if(!strcmp(word,"ifmore")) return cmdifmore(sim,cmd,line2);

	else if(!strcmp(word,"pointsource")) return cmdpointsource(sim,cmd,line2);
	else if(!strcmp(word,"volumesource")) return cmdvolumesource(sim,cmd,line2);
	else if(!strcmp(word,"killmol")) return cmdkillmol(sim,cmd,line2);
	else if(!strcmp(word,"killmolprob")) return cmdkillmolprob(sim,cmd,line2);
	else if(!strcmp(word,"killmolinsphere")) return cmdkillmolinsphere(sim,cmd,line2);
	else if(!strcmp(word,"equilmol")) return cmdequilmol(sim,cmd,line2);
	else if(!strcmp(word,"replacexyzmol")) return cmdreplacexyzmol(sim,cmd,line2);
	else if(!strcmp(word,"replacevolmol")) return cmdreplacevolmol(sim,cmd,line2);
	else if(!strcmp(word,"modulatemol")) return cmdmodulatemol(sim,cmd,line2);
	else if(!strcmp(word,"react1")) return cmdreact1(sim,cmd,line2);
	else if(!strcmp(word,"setrateint")) return cmdsetrateint(sim,cmd,line2);
	else if(!strcmp(word,"settimestep")) return cmdsettimestep(sim,cmd,line2);
	else if(!strcmp(word,"excludebox")) return cmdexcludebox(sim,cmd,line2);
	else if(!strcmp(word,"excludesphere")) return cmdexcludesphere(sim,cmd,line2);
	else if(!strcmp(word,"includeecoli")) return cmdincludeecoli(sim,cmd,line2);

	else if(!strcmp(word,"echo")) return cmdecho(sim,cmd,line2);
	else if(!strcmp(word,"molcount")) return cmdmolcount(sim,cmd,line2);
	else if(!strcmp(word,"listmols")) return cmdlistmols(sim,cmd,line2);
	else if(!strcmp(word,"listmols2")) return cmdlistmols2(sim,cmd,line2);
	else if(!strcmp(word,"listmols3")) return cmdlistmols3(sim,cmd,line2);
	else if(!strcmp(word,"molpos")) return cmdmolpos(sim,cmd,line2);
	else if(!strcmp(word,"molmoments")) return cmdmolmoments(sim,cmd,line2);
	else if(!strcmp(word,"savesim")) return cmdsavesim(sim,cmd,line2);
	else if(!strcmp(word,"meansqrdisp")) return cmdmeansqrdisp(sim,cmd,line2);
	SCMDCHECK(0,"command not recognized");
	return 1; }


/* ********** command routines *********** */

/****************** simulation control ********************/

int cmdstop(simptr sim,cmdptr cmd,char *line2) {
	return 2; }


int cmdpause(simptr sim,cmdptr cmd,char *line2) {
	char c;

	if(!sim->graphics) {
		fprintf(stderr,"Press enter: ");
		scanf("%c",&c);
		return 0; }
	gl2State(1);
	return 3; }


int cmdbeep(simptr sim,cmdptr cmd,char *line2) {
	fprintf(stderr,"\7");
	return 0; }


int cmdkeypress(simptr sim,cmdptr cmd,char *line2) {
	char c;
	int itct;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%c",&c);
	SCMDCHECK(itct==1,"cannot read character");
	gl2SetKeyPush((unsigned char) c);
	return 0; }


/****************** file manipulation ********************/


int cmdoverwrite(simptr sim,cmdptr cmd,char *line2) {
	SCMDCHECK(line2,"missing argument");
	SCMDCHECK(scmdoverwrite(sim->cmds,line2),"failed to open file");
	return docommand(sim,cmd,strnword(line2,2)); }


int cmdincrementfile(simptr sim,cmdptr cmd,char *line2) {
	SCMDCHECK(line2,"missing argument");
	SCMDCHECK(scmdincfile(sim->cmds,line2),"failed to open file");
	return docommand(sim,cmd,strnword(line2,2)); }



/****************** conditional ********************/


int cmdifno(simptr sim,cmdptr cmd,char *line2) {
	int itct,i,ll,flag,m;
	static char nm[STRCHAR];
	moleculeptr *mlist;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s",nm);
	SCMDCHECK(itct==1,"cannot read name");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	ll=(sim->mols->difc[i]==0);
	mlist=sim->mols->live[ll];
	flag=0;
	for(m=0;m<sim->mols->nl[ll]&&!flag;m++)
		if(mlist[m]->ident==i) flag++;
	if(!flag) return docommand(sim,cmd,strnword(line2,2));
	return 0; }


int cmdifless(simptr sim,cmdptr cmd,char *line2) {
	int itct,i,ll,ct,m,min;
	static char nm[STRCHAR];
	moleculeptr *mlist;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %i",nm,&min);
	SCMDCHECK(itct==2,"cannot read name");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	ll=(sim->mols->difc[i]==0);
	mlist=sim->mols->live[ll];
	ct=0;
	for(m=0;m<sim->mols->nl[ll]&&ct<min;m++)
		if(mlist[m]->ident==i) ct++;
	if(ct<min) return docommand(sim,cmd,strnword(line2,3));
	return 0; }


int cmdifmore(simptr sim,cmdptr cmd,char *line2) {
	int itct,i,ll,ct,m,max;
	static char nm[STRCHAR];
	moleculeptr *mlist;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %i",nm,&max);
	SCMDCHECK(itct==2,"cannot read name");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	ll=(sim->mols->difc[i]==0);
	mlist=sim->mols->live[ll];
	ct=0;
	for(m=0;m<sim->mols->nl[ll]&&ct<=max;m++)
		if(mlist[m]->ident==i) ct++;
	if(ct>max) return docommand(sim,cmd,strnword(line2,3));
	return 0; }


/****************** system manipulation ********************/


int cmdpointsource(simptr sim,cmdptr cmd,char *line2) {
	int itct,num,i,er,d;
	static char nm[STRCHAR];
	moleculeptr mptr;
	boxptr bptr;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %i",nm,&num);
	SCMDCHECK(itct==2,"read failure");
	SCMDCHECK(num>=0,"number cannot be negative");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	line2=strnword(line2,3);
	SCMDCHECK(line2,"missing location");
	itct=strreadnd(line2,sim->dim,sim->v1,NULL);
	SCMDCHECK(itct==sim->dim,"insufficient location dimensions");
	bptr=pos2box(sim,sim->v1);
	for(;num>0&&sim->mols->topd>0;num--) {
		mptr=sim->mols->dead[--sim->mols->topd];
		mptr->ident=i;
		for(d=0;d<sim->dim;d++) mptr->pos[d]=mptr->posx[d]=sim->v1[d];
		mptr->box=bptr; }
	er=molsort(sim->mols,0);
	if(er) return 2;
	SCMDCHECK(num==0,"not enough molecule spaces allocated");
	return 0; }


int cmdvolumesource(simptr sim,cmdptr cmd,char *line2) {
	int itct,num,i,er,dim,d;
	static char nm[STRCHAR];
	moleculeptr mptr;
	boxptr bptr;
	double *v1,*v2,*v3;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %i",nm,&num);
	SCMDCHECK(itct==2,"read failure");
	SCMDCHECK(num>=0,"number cannot be negative");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	line2=strnword(line2,3);
	SCMDCHECK(line2,"missing location");
	v1=sim->v1;
	v2=sim->v2;
	v3=sim->v3;
	dim=sim->dim;
	for(d=0;d<dim;d++) {
		SCMDCHECK(line2,"missing argument");
		itct=sscanf(line2,"%lf %lf",&v1[d],&v2[d]);
		SCMDCHECK(itct==2,"read failure");
		line2=strnword(line2,3); }

	for(;num>0&&sim->mols->topd>0;num--) {
		mptr=sim->mols->dead[--sim->mols->topd];
		mptr->ident=i;
		for(d=0;d<dim;d++) v3[d]=unirand(v1[d],v2[d]);
		bptr=pos2box(sim,v3);
		for(d=0;d<dim;d++) mptr->pos[d]=mptr->posx[d]=v3[d];
		mptr->box=bptr; }
	er=molsort(sim->mols,0);
	if(er) return 2;
	SCMDCHECK(num==0,"not enough molecule spaces allocated");
	return 0; }


int cmdkillmol(simptr sim,cmdptr cmd,char *line2) {
	int itct,i,ll,m,er;
	static char nm[STRCHAR];
	moleculeptr *mlist;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s",nm);
	SCMDCHECK(itct==1,"cannot read name");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	ll=(sim->mols->difc[i]==0);
	mlist=sim->mols->live[ll];
	for(m=0;m<sim->mols->nl[ll];m++)
		if(mlist[m]->ident==i) mlist[m]->ident=0;
	er=molsort(sim->mols,0);
	if(er) return 2;
	return 0; }


int cmdkillmolprob(simptr sim,cmdptr cmd,char *line2) {
	int itct,i,ll,m,er;
	static char nm[STRCHAR];
	moleculeptr *mlist;
	double prob;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %lf",nm,&prob);
	SCMDCHECK(itct==2,"killmolprob format: name probability");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	SCMDCHECK(prob>=0&&prob<=1,"probability needs to be between 0 and 1");
	ll=(sim->mols->difc[i]==0);
	mlist=sim->mols->live[ll];
	for(m=0;m<sim->mols->nl[ll];m++)
		if(mlist[m]->ident==i&&coinrand(prob)) mlist[m]->ident=0;
	er=molsort(sim->mols,0);
	if(er) return 2;
	return 0; }


int cmdkillmolinsphere(simptr sim,cmdptr cmd,char *line2) {
	int itct,i,s,ll,m,er;
	static char nm1[STRCHAR],nm2[STRCHAR];
	moleculeptr *mlist;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %s",nm1,nm2);
	SCMDCHECK(itct==2,"cannot read molecule name or surface name");
	if(!strcmp(nm1,"all")) {
		i=-1;
		ll=0;
		mlist=NULL; }
	else {
		i=stringfind(sim->name,sim->nident,nm1);
		SCMDCHECK(i>=1,"name not recognized");
		ll=(sim->mols->difc[i]==0);
		mlist=sim->mols->live[ll]; }

	if(!strcmp(nm2,"all")) s=-1;
	else {
		if(!sim->srfss) return 0;
		s=stringfind(sim->srfss->snames,sim->srfss->nsrf,nm2);
		SCMDCHECK(s>=0,"surface not recognized"); }

	if(i>0) {				// certain identity
		for(m=0;m<sim->mols->nl[ll];m++)
			if(mlist[m]->ident==i)
				if(molinpanels(sim,ll,m,s,'s')) mlist[m]->ident=0; }
	else {					// all identities
		for(ll=0;ll<2;ll++) {
			mlist=sim->mols->live[ll];
			for(m=0;m<sim->mols->nl[ll];m++)
				if(molinpanels(sim,ll,m,s,'s')) mlist[m]->ident=0; }}

	er=molsort(sim->mols,0);
	if(er) return 2;
	return 0; }


int cmdequilmol(simptr sim,cmdptr cmd,char *line2) {
	int itct,i1,i2,m,er,ll1,ll2;
	static char nm[STRCHAR],nm2[STRCHAR];
	moleculeptr *mlist;
	double prob;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %s %lf",nm,nm2,&prob);
	SCMDCHECK(itct==3,"read failure");
	SCMDCHECK(prob>=0&&prob<=1,"probability is out of bounds");
	i1=stringfind(sim->name,sim->nident,nm);
	i2=stringfind(sim->name,sim->nident,nm2);
	SCMDCHECK(i1>=1&&i2>=1,"name not recognized");
	ll1=sim->mols->difc[i1]==0;
	ll2=sim->mols->difc[i2]==0;
	if(ll1||ll2) {
		mlist=sim->mols->live[1];
		for(m=0;m<sim->mols->nl[1];m++)
			if(mlist[m]->ident==i1||mlist[m]->ident==i2)
				mlist[m]->ident=coinrand(prob)?i2:i1; }
	if(!ll1||!ll2) {
		mlist=sim->mols->live[0];
		for(m=0;m<sim->mols->nl[0];m++)
			if(mlist[m]->ident==i1||mlist[m]->ident==i2)
				mlist[m]->ident=coinrand(prob)?i2:i1; }
	if(ll1!=ll2) {
		er=molsort(sim->mols,1);
		if(er) return 2; }
	return 0; }


int cmdreplacexyzmol(simptr sim,cmdptr cmd,char *line2) {
	int itct,i1,m,d,dim,er;
	static char nm[STRCHAR];
	moleculeptr *mlist;
	boxptr bptr;
	double *v1;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s",nm);
	SCMDCHECK(itct==1,"cannot read name");
	i1=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i1>=1,"name not recognized");
	dim=sim->dim;
	v1=sim->v1;
	line2=strnword(line2,2);
	itct=strreadnd(line2,dim,v1,NULL);
	SCMDCHECK(itct==dim,"insufficient dimensions entered");
	bptr=pos2box(sim,v1);
	mlist=bptr->mol[1];
	for(m=0;m<bptr->nmol[1];m++) {
		for(d=0;d<dim;d++)
			if(mlist[m]->pos[d]!=v1[d]) d=dim+1;
		if(d==dim) {
			mlist[m]->ident=i1;
			m=bptr->nmol[1]+1; }}
	if(m>bptr->nmol[1]&&sim->mols->difc[i1]!=0) {
		er=molsort(sim->mols,1);
		if(er) return 2; }
	return 0; }


int cmdreplacevolmol(simptr sim,cmdptr cmd,char *line2) {
	int m,itct,dim,d,b,b1,b2,i1,i2,ll;
	static char nm1[STRCHAR],nm2[STRCHAR];
	double *pos,*v1,*v2,frac;
	boxptr bptr1,bptr2,bptr;
	boxssptr boxs;
	moleculeptr *mlist;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %s %lf",nm1,nm2,&frac);
	SCMDCHECK(itct==3,"cannot read names");
	i1=stringfind(sim->name,sim->nident,nm1);
	SCMDCHECK(i1>=1,"first name not recognized");
	i2=stringfind(sim->name,sim->nident,nm2);
	SCMDCHECK(i2>=1,"second name not recognized");
	SCMDCHECK(frac>=0&&frac<=1,"fraction out of bounds");
	line2=strnword(line2,4);
	dim=sim->dim;
	boxs=sim->boxs;
	v1=sim->v1;
	v2=sim->v2;
	for(d=0;d<dim;d++) {
		SCMDCHECK(line2,"missing argument");
		itct=sscanf(line2,"%lf %lf",&v1[d],&v2[d]);
		SCMDCHECK(itct==2,"read failure");
		line2=strnword(line2,3); }
	bptr1=pos2box(sim,v1);
	bptr2=pos2box(sim,v2);
	SCMDCHECK(bptr1&&bptr2,"failed to get pointer to virtual box");
	b1=indx2addZV(bptr1->indx,boxs->side,dim);
	b2=indx2addZV(bptr2->indx,boxs->side,dim);
	ll=(sim->mols->difc[i1]==0);
	for(b=b1;b<=b2;b=nextaddZV(b,bptr1->indx,bptr2->indx,boxs->side,dim)) {
		bptr=boxs->blist[b];
		mlist=bptr->mol[ll];
		for(m=0;m<bptr->nmol[ll];m++) {
			if(mlist[m]->ident==i1) {
				pos=mlist[m]->pos;
				for(d=0;d<dim;d++) if(pos[d]<v1[d]||pos[d]>v2[d]) d=dim+1;
				if(d==dim&&coinrand(frac)) mlist[m]->ident=i2; }}}
	return 0; }


int cmdmodulatemol(simptr sim,cmdptr cmd,char *line2) {
	int itct,i1,i2,m,er,ll1,ll2;
	static char nm[STRCHAR],nm2[STRCHAR];
	moleculeptr *mlist;
	double freq,shift,prob;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %s %lf %lf",nm,nm2,&freq,&shift);
	SCMDCHECK(itct==4,"read failure");
	i1=stringfind(sim->name,sim->nident,nm);
	i2=stringfind(sim->name,sim->nident,nm2);
	SCMDCHECK(i1>=1&&i2>=1,"name not recognized");
	ll1=sim->mols->difc[i1]==0;
	ll2=sim->mols->difc[i2]==0;
	prob=0.5*(1.0-cos(freq*sim->time+shift));
	if(ll1||ll2) {
		mlist=sim->mols->live[1];
		for(m=0;m<sim->mols->nl[1];m++)
			if(mlist[m]->ident==i1||mlist[m]->ident==i2)
				mlist[m]->ident=coinrand(prob)?i2:i1; }
	if(!ll1||!ll2) {
		mlist=sim->mols->live[0];
		for(m=0;m<sim->mols->nl[0];m++)
			if(mlist[m]->ident==i1||mlist[m]->ident==i2)
				mlist[m]->ident=coinrand(prob)?i2:i1; }
	if(ll1!=ll2) {
		er=molsort(sim->mols,1);
		if(er) return 2; }
	return 0; }


int cmdreact1(simptr sim,cmdptr cmd,char *line2) {
	int itct,i,ll,m,r,er;
	static char nm[STRCHAR],rnm[STRCHAR];
	moleculeptr *mlist;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %s",nm,rnm);
	SCMDCHECK(itct==2,"read failure");
	SCMDCHECK(sim->rxn[1],"no first order reactions defined");
	i=stringfind(sim->name,sim->nident,nm);
	r=stringfind(sim->rxn[1]->rname,sim->rxn[1]->total,rnm);
	SCMDCHECK(r>=0,"reaction not recognized");
	SCMDCHECK(i>=1,"molecule not recognized");
	ll=(sim->mols->difc[i]==0);
	mlist=sim->mols->live[ll];
	for(m=0;m<sim->mols->nl[ll];m++)
		if(mlist[m]->ident==i)
			if(doreact(sim->rxn[1],r,mlist[m],NULL,sim)) return 1;
	er=molsort(sim->mols,0);
	if(er) return 2;
	return 0; }


int cmdsetrateint(simptr sim,cmdptr cmd,char *line2) {
	int itct,r,order;
	static char rnm[STRCHAR];
	double rateint;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %lf",rnm,&rateint);
	SCMDCHECK(itct==2,"read failure");
	r=-1;
	if(sim->rxn[0]) r=stringfind(sim->rxn[0]->rname,sim->rxn[0]->total,rnm);
	if(r>=0) order=0;
	else {
		if(sim->rxn[1]) r=stringfind(sim->rxn[1]->rname,sim->rxn[1]->total,rnm);
		if(r>=0) order=1;
		else {
			if(sim->rxn[2]) r=stringfind(sim->rxn[2]->rname,sim->rxn[2]->total,rnm);
			if(r>=0) order=2;
			else SCMDCHECK(0,"reaction not recognized"); }}
	SCMDCHECK(rateint>=0,"internal rate cannot be negative");
	if(order<2) sim->rxn[order]->rate2[r]=rateint;
	else sim->rxn[order]->rate2[r]=rateint*rateint;
	return 0; }


int cmdsettimestep(simptr sim,cmdptr cmd,char *line2) {
	int itct;
	double dt;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%lf",&dt);
	SCMDCHECK(itct==1,"read failure");
	SCMDCHECK(dt>0,"time step must be >0");
	sim->dt=dt;
	setdiffusion(sim);
	return 0; }


int cmdexcludebox(simptr sim,cmdptr cmd,char *line2) {
	int m,itct,dim,d,b,b1,b2;
	double *pos,*v1,*v2;
	boxptr bptr1,bptr2,bptr;
	boxssptr boxs;
	moleculeptr *mlist;

	dim=sim->dim;
	boxs=sim->boxs;
	v1=sim->v1;
	v2=sim->v2;
	for(d=0;d<dim;d++) {
		SCMDCHECK(line2,"missing argument");
		itct=sscanf(line2,"%lf %lf",&v1[d],&v2[d]);
		SCMDCHECK(itct==2,"read failure");
		line2=strnword(line2,3); }
	bptr1=pos2box(sim,v1);
	bptr2=pos2box(sim,v2);
	SCMDCHECK(bptr1&&bptr2,"failed to get pointer to virtual box");
	b1=indx2addZV(bptr1->indx,boxs->side,dim);
	b2=indx2addZV(bptr2->indx,boxs->side,dim);
	for(b=b1;b<=b2;b=nextaddZV(b,bptr1->indx,bptr2->indx,boxs->side,dim)) {
		bptr=boxs->blist[b];
		mlist=bptr->mol[0];
		for(m=0;m<bptr->nmol[0];m++) {
			pos=mlist[m]->pos;
			for(d=0;d<dim;d++) if(pos[d]<v1[d]||pos[d]>v2[d]) d=dim+1;
			if(d==dim) {
				pos=mlist[m]->posx;
				for(d=0;d<dim;d++) if(pos[d]<v1[d]||pos[d]>v2[d]) d=dim+1;
				if(d>dim) copyVD(mlist[m]->posx,mlist[m]->pos,dim); }}}
	return 0; }


int cmdexcludesphere(simptr sim,cmdptr cmd,char *line2) {
	int m,itct,dim,d,b,b1,b2;
	double *pos,*v1,*v2,*v3,rad,dist;
	boxptr bptr1,bptr2,bptr;
	boxssptr boxs;
	moleculeptr *mlist;

	dim=sim->dim;
	boxs=sim->boxs;
	v1=sim->v1;
	v2=sim->v2;
	v3=sim->v3;
	for(d=0;d<dim;d++) {
		SCMDCHECK(line2,"missing argument");
		itct=sscanf(line2,"%lf",&v3[d]);
		SCMDCHECK(itct==1,"read failure");
		line2=strnword(line2,2); }
	SCMDCHECK(line2,"missing radius");
	itct=sscanf(line2,"%lf",&rad);
	SCMDCHECK(itct==1,"read failure");
	dist=rad*sqrt(dim);
	for(d=0;d<dim;d++) {
		v1[d]=v3[d]-dist;
		v2[d]=v3[d]+dist; }
	rad*=rad;
	bptr1=pos2box(sim,v1);
	bptr2=pos2box(sim,v2);
	SCMDCHECK(bptr1&&bptr2,"failed to get pointer to virtual box");
	b1=indx2addZV(bptr1->indx,boxs->side,dim);
	b2=indx2addZV(bptr2->indx,boxs->side,dim);
	for(b=b1;b<=b2;b=nextaddZV(b,bptr1->indx,bptr2->indx,boxs->side,dim)) {
		bptr=boxs->blist[b];
		mlist=bptr->mol[0];
		for(m=0;m<bptr->nmol[0];m++) {
			pos=mlist[m]->pos;
			for(dist=0,d=0;d<dim;d++) if((dist+=(pos[d]-v3[d])*(pos[d]-v3[d]))>rad) d=dim+1;
			if(d==dim) {
				pos=mlist[m]->posx;
				for(dist=0,d=0;d<dim;d++) if((dist+=(pos[d]-v3[d])*(pos[d]-v3[d]))>rad) d=dim+1;
				if(d>dim) copyVD(mlist[m]->posx,mlist[m]->pos,dim); }}}
	return 0; }


int cmdincludeecoli(simptr sim,cmdptr cmd,char *line2) {
	int m,nl;
	moleculeptr *mlist;
	double rad,length,*v1;
	wallptr *wlist;

	SCMDCHECK(sim->dim==3,"system is not 3 dimensional");
	v1=sim->v1;
	wlist=sim->wlist;
	rad=0.5*(wlist[3]->pos-wlist[2]->pos);
	length=wlist[1]->pos-wlist[0]->pos;
	v1[0]=wlist[0]->pos;
	v1[1]=0.5*(wlist[2]->pos+wlist[3]->pos);
	v1[2]=0.5*(wlist[4]->pos+wlist[5]->pos);
	mlist=sim->mols->live[0];
	nl=sim->mols->nl[0];
	for(m=0;m<nl;m++)
		if(!insideecoli(mlist[m]->pos,v1,rad,length)) {
			if(insideecoli(mlist[m]->posx,v1,rad,length)) copyVD(mlist[m]->posx,mlist[m]->pos,3);
			else putinecoli(mlist[m]->pos,v1,rad,length); }
	return 0; }


/**************** observation commands ******************/


int cmdecho(simptr sim,cmdptr cmd,char *line2) {
	FILE *fptr;
	char *termqt,str[STRCHAR];

	fptr=scmdgetfptr(sim->cmds,line2);
	SCMDCHECK(fptr,"file name not recognized");
	line2=strnword(line2,2);
	SCMDCHECK(line2=strchr(line2,'"'),"no starting quote on string");
	strncpy(str,line2+1,STRCHAR);
	SCMDCHECK(termqt=strchr(str,'"'),"no terminal quote on string");
	*termqt='\0';
	strbslash2escseq(str);
	fprintf(fptr,"%s",str);
	return 0; }


int cmdmolcount(simptr sim,cmdptr cmd,char *line2) {
	FILE *fptr;
	int ll,m,*ct;

	fptr=scmdgetfptr(sim->cmds,line2);
	SCMDCHECK(fptr,"file name not recognized");
	ct=sim->z1;
	setstdZV(ct,sim->nident,0);
	for(ll=0;ll<2;ll++)
		for(m=0;m<sim->mols->nl[ll];m++)
			ct[sim->mols->live[ll][m]->ident]++;
	fprintf(fptr,"%1.5f ",sim->time);
	fprintZV(fptr,ct+1,sim->nident-1);
	return 0; }


int cmdlistmols(simptr sim,cmdptr cmd,char *line2) {
	int m,ll;
	moleculeptr mptr;
	FILE *fptr;

	fptr=scmdgetfptr(sim->cmds,line2);
	SCMDCHECK(fptr,"file name not recognized");
	for(ll=0;ll<2;ll++)
		for(m=0;m<sim->mols->nl[ll];m++) {
			mptr=sim->mols->live[ll][m];
			fprintf(fptr,"%s ",sim->name[mptr->ident]);
			fprintVD(fptr,mptr->pos,sim->dim); }
	return 0; }


int cmdlistmols2(simptr sim,cmdptr cmd,char *line2) {
	int m,ll,invk;
	moleculeptr mptr;
	FILE *fptr;

	fptr=scmdgetfptr(sim->cmds,line2);
	SCMDCHECK(fptr,"file name not recognized");
	invk=cmd?cmd->invoke:0;
	for(ll=0;ll<2;ll++)
		for(m=0;m<sim->mols->nl[ll];m++) {
			mptr=sim->mols->live[ll][m];
			fprintf(fptr,"%i %i ",invk,mptr->ident);
			fprintVD(fptr,mptr->pos,sim->dim); }
	return 0; }


int cmdlistmols3(simptr sim,cmdptr cmd,char *line2) {
	int i,itct,m,ll,dim,invk;
	moleculeptr *mlist,mptr;
	FILE *fptr;
	static char nm[STRCHAR];

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s",nm);
	SCMDCHECK(itct==1,"cannot read name");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	line2=strnword(line2,2);
	fptr=scmdgetfptr(sim->cmds,line2);
	SCMDCHECK(fptr,"file name not recognized");
	invk=cmd?cmd->invoke:0;

	dim=sim->dim;
	ll=(sim->mols->difc[i]==0);
	mlist=sim->mols->live[ll];
	for(m=0;m<sim->mols->nl[ll];m++)
		if(mlist[m]->ident==i) {
			mptr=mlist[m];
			fprintf(fptr,"%i %i ",invk,mptr->ident);
			fprintVD(fptr,mptr->pos,sim->dim); }
	return 0; }


int cmdmolpos(simptr sim,cmdptr cmd,char *line2) {
	int i,d,itct,m,ll,dim;
	moleculeptr *mlist;
	FILE *fptr;
	static char nm[STRCHAR];

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s",nm);
	SCMDCHECK(itct==1,"cannot read name");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	line2=strnword(line2,2);
	fptr=scmdgetfptr(sim->cmds,line2);
	SCMDCHECK(fptr,"file name not recognized");

	dim=sim->dim;
	ll=(sim->mols->difc[i]==0);
	mlist=sim->mols->live[ll];
	fprintf(fptr,"%lf ",sim->time);
	for(m=0;m<sim->mols->nl[ll];m++)
		if(mlist[m]->ident==i)
			for(d=0;d<dim;d++)
				fprintf(fptr,"%lf ",mlist[m]->pos[d]);
	fprintf(fptr,"\n");
	return 0; }


int cmdmolmoments(simptr sim,cmdptr cmd,char *line2) {
	static char nm[STRCHAR];
	int i,itct,ll,m,ctr,dim,d,d2;
	double *v1,*m1;
	FILE *fptr;
	moleculeptr *mlist;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s",nm);
	SCMDCHECK(itct==1,"cannot read name");
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	line2=strnword(line2,2);
	fptr=scmdgetfptr(sim->cmds,line2);
	SCMDCHECK(fptr,"file name not recognized");

	dim=sim->dim;
	ll=(sim->mols->difc[i]==0);
	mlist=sim->mols->live[ll];
	v1=sim->v1;
	m1=sim->m1;
	ctr=0;
	for(d=0;d<dim;d++) v1[d]=0;
	for(d=0;d<dim*dim;d++) m1[d]=0;
	for(m=0;m<sim->mols->nl[ll];m++)
		if(mlist[m]->ident==i) {
			ctr++;
			for(d=0;d<dim;d++) v1[d]+=mlist[m]->pos[d]; }
	for(d=0;d<dim;d++) v1[d]/=1.0*ctr;
	for(m=0;m<sim->mols->nl[ll];m++)
		if(mlist[m]->ident==i) {
			for(d=0;d<dim;d++)
				for(d2=0;d2<dim;d2++)
					m1[d*dim+d2]+=(mlist[m]->pos[d]-v1[d])*(mlist[m]->pos[d2]-v1[d2]); }
	fprintf(fptr,"%lf %i",sim->time,ctr);
	for(d=0;d<dim;d++) fprintf(fptr," %lf",v1[d]);
	for(d=0;d<dim;d++)
		for(d2=0;d2<dim;d2++)
			fprintf(fptr," %lf",m1[d*dim+d2]/ctr);
	fprintf(fptr,"\n");
	return 0; }


int cmdsavesim(simptr sim,cmdptr cmd,char *line2) {
	int i,i1,j,nident,dim,d,ll,m,r,p,order,s;
	FILE *fptr;
	char **name,**rname;
	cmdptr cmd2;
	void *voidptr;
	rxnptr rxn;
	surfaceptr srf;
	panelptr pnl;

	fptr=scmdgetfptr(sim->cmds,line2);
	SCMDCHECK(fptr,"file name not recognized");
	dim=sim->dim;
	nident=sim->nident;
	name=sim->name;
	fprintf(fptr,"# Configuration file automatically created by Smoldyn\n\n");
	if(sim->graphics==0) fprintf(fptr,"graphics none\n");
	else if(sim->graphics==1) fprintf(fptr,"graphics opengl\n");
	else if(sim->graphics==2) fprintf(fptr,"graphics opengl_good\n");
	fprintf(fptr,"graphic_iter %i\n",sim->graphicit);
	if(sim->tiffit>0) {
		fprintf(fptr,"tiff_iter %i\n",sim->tiffit);
		fprintf(fptr,"tiff_name %s\n",gl2SetOptionStr("TiffName",NULL));
		fprintf(fptr,"tiff_min %i\n",gl2SetOptionInt("TiffNumber",-1));
		fprintf(fptr,"tiff_max %i\n",gl2SetOptionInt("TiffNumMax",-1)); }
	fprintf(fptr,"\n");
	fprintf(fptr,"dim %i\n",dim);
	fprintf(fptr,"max_names %i\n",sim->maxident);
	for(i=1;i<nident;i++) fprintf(fptr,"name %s\n",name[i]);
	fprintf(fptr,"max_mol %i\n",sim->mols->max);
	if(sim->boxs->mpbox) fprintf(fptr,"molperbox %lf\n",sim->boxs->mpbox);
	else if(sim->boxs->boxsize) fprintf(fptr,"boxsize %lf\n",sim->boxs->boxsize);
	fprintf(fptr,"accuracy %lf\n",sim->accur);
	fprintf(fptr,"frame_thickness %lf\n",sim->framepts);
	fprintf(fptr,"frame_color %lf %lf %lf %lf\n",sim->framecolor[0],sim->framecolor[1],sim->framecolor[2],sim->framecolor[3]);
	fprintf(fptr,"background_color %lf %lf %lf %lf\n",sim->backcolor[0],sim->backcolor[1],sim->backcolor[2],sim->backcolor[3]);

	for(i=1;i<nident;i++)
		fprintf(fptr,"difc %s %lf\n",name[i],sim->mols->difc[i]);
	for(i=1;i<nident;i++)
		if(sim->mols->difm[i]) {
			fprintf(fptr,"difm %s",name[i]);
			for(d=0;d<dim*dim;d++)
				fprintf(fptr," %lf",sim->mols->difm[i][d]);
			fprintf(fptr,"\n"); }
	for(i=1;i<nident;i++)
		fprintf(fptr,"color %s %lf %lf %lf\n",name[i],sim->mols->color[i][0],sim->mols->color[i][1],sim->mols->color[i][2]);
	for(i=1;i<nident;i++)
		fprintf(fptr,"display_size %s %lf\n",name[i],sim->mols->display[i]);
	fprintf(fptr,"\ntime_start %lf\n",sim->tmin);
	fprintf(fptr,"time_stop %lf\n",sim->tmax);
	fprintf(fptr,"time_step %lf\n",sim->dt);
	fprintf(fptr,"time_now %lf\n",sim->time);
	fprintf(fptr,"# rand_seed for prior simulation was %u\n",sim->randseed);
	fprintf(fptr,"rand_seed %i  # this is a new random number\n\n",rand());
	for(d=0;d<dim;d++) {
		fprintf(fptr,"low_wall %i %lf %c\n",d,sim->wlist[2*d]->pos,sim->wlist[2*d]->type);
		fprintf(fptr,"high_wall %i %lf %c\n",d,sim->wlist[2*d+1]->pos,sim->wlist[2*d+1]->type); }
	fprintf(fptr,"\n");

	if(sim->srfss) {
		fprintf(fptr,"max_surface %i\n",sim->srfss->maxsrf);
		for(s=0;s<sim->srfss->nsrf;s++) {
			srf=sim->srfss->srflist[s];
			fprintf(fptr,"start_surface\n");
			fprintf(fptr,"name %s\n",sim->srfss->snames[s]);
			for(i=1;i<nident;i++) {
				fprintf(fptr,"action_front %s %c\n",name[i],srf->faction[i]);
				fprintf(fptr,"action_back %s %c\n",name[i],srf->baction[i]); }
			fprintf(fptr,"color_front %lf %lf %lf %lf\n",srf->fcolor[0],srf->fcolor[1],srf->fcolor[2],srf->fcolor[3]);
			fprintf(fptr,"color_back %lf %lf %lf %lf\n",srf->bcolor[0],srf->bcolor[1],srf->bcolor[2],srf->bcolor[3]);
			fprintf(fptr,"thickness %lf\n",srf->edgepts);
			fprintf(fptr,"polygon_front %c\n",srf->fpolymode);
			fprintf(fptr,"polygon_back %c\n",srf->bpolymode);
			if(srf->maxpanel[0]) fprintf(fptr,"max_panels r %i\n",srf->maxpanel[0]);
			if(srf->maxpanel[1]) fprintf(fptr,"max_panels t %i\n",srf->maxpanel[1]);
			if(srf->maxpanel[2]) fprintf(fptr,"max_panels s %i\n",srf->maxpanel[2]);
			if(srf->maxpanel[3]) fprintf(fptr,"max_panels c %i\n",srf->maxpanel[3]);
			if(srf->maxpanel[4]) fprintf(fptr,"max_panels h %i\n",srf->maxpanel[4]);
			for(p=0;p<srf->npanel[0];p++) {
				pnl=srf->panels[0][p];
				fprintf(fptr,"panel r %c%i",pnl->front[0]>0?'+':'-',(int)(pnl->front[1]));
				for(d=0;d<dim;d++) fprintf(fptr," %lf",pnl->point[0][d]);
				for(d=0;d<dim;d++)
					if(d!=pnl->front[1]) fprintf(fptr," %lf",pnl->point[dim==2?1:2][d]-pnl->point[0][d]);
				fprintf(fptr,"\n"); }
			for(p=0;p<srf->npanel[1];p++) {
				pnl=srf->panels[1][p];
				fprintf(fptr,"panel t");
				for(j=0;j<dim;j++)
					for(d=0;d<dim;d++)
						fprintf(fptr," %lf",pnl->point[j][d]);
				fprintf(fptr,"\n"); }
			for(p=0;p<srf->npanel[2];p++) {
				pnl=srf->panels[2][p];
				fprintf(fptr,"panel s");
				for(d=0;d<dim;d++) fprintf(fptr," %lf",pnl->point[0][d]);
				fprintf(fptr," %c%lf",pnl->front[0]>0?'+':'-',pnl->point[1][0]);
				if(dim>1) fprintf(fptr," %lf",pnl->point[1][1]);
				if(dim>2) fprintf(fptr," %lf",pnl->point[1][2]);
				fprintf(fptr,"\n"); }
			for(p=0;p<srf->npanel[3];p++) {
				pnl=srf->panels[3][p];
				fprintf(fptr,"panel c");
				for(d=0;d<dim;d++) fprintf(fptr," %lf",pnl->point[0][d]);
				for(d=0;d<dim;d++) fprintf(fptr," %lf",pnl->point[1][d]);
				fprintf(fptr," %c%lf",pnl->front[dim==2?2:0]>0?'+':'-',pnl->point[2][0]);
				if(dim>2) fprintf(fptr," %lf %lf",pnl->point[2][1],pnl->point[2][2]);
				fprintf(fptr,"\n"); }
			for(p=0;p<srf->npanel[4];p++) {
				pnl=srf->panels[4][p];
				fprintf(fptr,"panel h");
				for(d=0;d<dim;d++) fprintf(fptr," %lf",pnl->point[0][d]);
				fprintf(fptr," %c%lf",pnl->front[0]>0?'+':'-',pnl->point[1][0]);
				for(d=0;d<dim;d++) fprintf(fptr," %lf",pnl->point[2][d]);
				if(dim>=2) fprintf(fptr," %lf",pnl->point[1][1]);
				if(dim>2) fprintf(fptr," %lf",pnl->point[1][2]);
				fprintf(fptr,"\n"); }
			fprintf(fptr,"end_surface\n\n"); }}

	if(strlen(sim->cmds->froot)||sim->cmds->nfile) fprintf(fptr,"\n");
	if(strlen(sim->cmds->froot))
		fprintf(fptr,"output_root %s\n",sim->cmds->froot);
	if(sim->cmds->nfile) {
		fprintf(fptr,"output_files");
		for(i=0;i<sim->cmds->nfile;i++)
			fprintf(fptr," %s",sim->cmds->fname[i]);
		fprintf(fptr,"\n"); }
	for(i=0;i<sim->cmds->nfile;i++)
		if(sim->cmds->fsuffix[i])
			fprintf(fptr,"output_file_number %s %i\n",sim->cmds->fname[i],sim->cmds->fsuffix[i]);
	fprintf(fptr,"\nmax_cmd %i\n",qmaxlength(sim->cmds->cmd)-1);
	i=-1;
	voidptr=NULL;
	while((i=qnext(i,NULL,&voidptr,sim->cmds->cmd))>=0) {
		cmd2=(cmdptr)voidptr;
		fprintf(fptr,"cmd i %lf %lf %lf %s\n",cmd2->on,cmd2->off,cmd2->dt,cmd2->str); }
	fprintf(fptr,"\n");

	for(order=0;order<=2;order++)
		if(sim->rxn[order]) {
			rxn=sim->rxn[order];
			rname=rxn->rname;
			fprintf(fptr,"start_reaction\n");
			fprintf(fptr,"order %i\n",order);
			fprintf(fptr,"max_rxn %i\n",rxn->total);
			if(order==1) {
				for(i=0;i<nident;i++)
					if(rxn->nrxn[i]) {
						fprintf(fptr,"reactant %s %s",name[i],rname[rxn->table[i][0]]);
						for(j=1;j<rxn->nrxn[i];j++) fprintf(fptr," + %s",rname[rxn->table[i][j]]);
						fprintf(fptr,"\n"); }}
			else if(order==2)
				for(i=0;i<nident;i++)
					for(i1=0;i1<=i;i1++)
						if(rxn->nrxn[i*nident+i1]) {
							fprintf(fptr,"reactant %s + %s %s",name[i],name[i1],rname[rxn->table[i*nident+i1][0]]);
							for(j=1;j<rxn->nrxn[i*nident+i1];j++) fprintf(fptr," + %s",rname[rxn->table[i*nident+i1][j]]);
							fprintf(fptr,"\n"); }
			for(r=0;r<rxn->total;r++) {
				if(rxn->rate[r]>=0) fprintf(fptr,"rate %s %lf\n",rname[r],rxn->rate[r]);
				else if(rxn->rate2[r]>=0) fprintf(fptr,"rate_internal %s %lf\n",rname[r],order<2?rxn->rate2[r]:sqrt(rxn->rate2[r]));
				if(rxn->nprod[r])
					fprintf(fptr,"product %s %s",rname[r],name[rxn->prod[r][0]->ident]);
				for(p=1;p<rxn->nprod[r];p++)
					fprintf(fptr," + %s",name[rxn->prod[r][p]->ident]);
				fprintf(fptr,"\n");
				if(rxn->rpart[r]=='i') fprintf(fptr,"product_param %s i\n",rname[r]);
				else if(strchr("pxrbqys",rxn->rpart[r])) fprintf(fptr,"product_param %s %c %lf\n",rname[r],rxn->rpart[r],rxn->rpar[r]);
				else if(strchr("of",rxn->rpart[r])) {
					for(p=0;p<rxn->nprod[r];p++) {
						fprintf(fptr,"product_param %s %c %s\n",rname[r],rxn->rpart[r],name[rxn->prod[r][p]->ident]);
						for(d=0;d<dim;d++) fprintf(fptr," %lf",rxn->prod[r][p]->pos[d]);
						fprintf(fptr,"\n"); }}}
			fprintf(fptr,"end_reaction\n\n"); }
	for(ll=0;ll<2;ll++)
		for(m=0;m<sim->mols->nl[ll];m++) {
			fprintf(fptr,"mol 1 %s ",name[sim->mols->live[ll][m]->ident]);
			for(d=0;d<dim;d++) fprintf(fptr,"%lf ",sim->mols->live[ll][m]->pos[d]);
			fprintf(fptr,"\n"); }
	fprintf(fptr,"\nend_file\n");
	return 0; }


void cmdmeansqrdispfree(cmdptr cmd) {
	int j;

	if(cmd->v3&&cmd->i1)
		for(j=0;j<cmd->i1;j++)
			free(((int**)(cmd->v3))[j]);
	if(cmd->v3) free(cmd->v3);
	if(cmd->v2&&cmd->i1)
		for(j=0;j<cmd->i1;j++)
			free(((double**)(cmd->v2))[j]);
	if(cmd->v2) free(cmd->v2);
	if(cmd->v1) free(cmd->v1);
	return; }


int cmdmeansqrdisp(simptr sim,cmdptr cmd,char *line2) {
	static char nm[STRCHAR],dimstr[STRCHAR];
	static double space[5];
	int i,j,d,itct,ll,dim,ctr,m,msddim;
	FILE *fptr;
	moleculeptr *mlist;
	double sum,diff,*pos,**v2;
	long int *v1;

	SCMDCHECK(line2,"missing argument");
	itct=sscanf(line2,"%s %s",nm,dimstr);
	SCMDCHECK(itct==2,"cannot read name and/or dimension");
	if(!strcmp(dimstr,"all")) msddim=-1;
	else {
		itct=sscanf(dimstr,"%i",&msddim);
		SCMDCHECK(itct==1,"cannot read dimension");
		SCMDCHECK(msddim>=0&&msddim<sim->dim,"dimension out of range"); }
	i=stringfind(sim->name,sim->nident,nm);
	SCMDCHECK(i>=1,"name not recognized");
	line2=strnword(line2,3);
	fptr=scmdgetfptr(sim->cmds,line2);
	SCMDCHECK(fptr,"file name not recognized");

	SCMDCHECK(cmd->i2!=2,"error on setup");					// failed before, don't try again

	ll=sim->mols->difc[i]==0;
	dim=sim->dim;
	mlist=sim->mols->live[ll];
	SCMDCHECK(dim<=5,"memory is allocated for more than 5 dimensional space");

	if(!cmd->i2) {										// test for first run
		cmd->i2=1;											// if first run, set up data structures
		for(d=0;d<dim;d++)
			space[d]=sim->wlist[2*d+1]->pos-sim->wlist[2*d]->pos;
		ctr=0;
		for(m=0;m<sim->mols->nl[ll];m++)
			if(mlist[m]->ident==i) ctr++;
		cmd->i1=ctr;										// size of arrays
		cmd->freefn=&cmdmeansqrdispfree;
		cmd->v1=calloc(ctr,sizeof(long int));	// v1 is serial numbers
		if(!cmd->v1) {cmd->i2=2;return 1;}
		for(j=0;j<ctr;j++) ((long int*)cmd->v1)[j]=0;
		cmd->v2=calloc(ctr,sizeof(double**));	// v2 is positions
		if(!cmd->v2) {cmd->i2=2;return 1;}
		for(j=0;j<ctr;j++) ((double**)cmd->v2)[j]=NULL;
		for(j=0;j<ctr;j++) {
			((double**)cmd->v2)[j]=(double*)calloc(dim,sizeof(double));
			if(!((double**)cmd->v2)[j]) {cmd->i2=2;return 1;}
			for(d=0;d<dim;d++) ((double**)cmd->v2)[j][d]=0; }
		ctr=0;
		for(m=0;m<sim->mols->nl[ll];m++) {
			if(mlist[m]->ident==i) {
				((long int*)cmd->v1)[ctr]=mlist[m]->serno;
				for(d=0;d<dim;d++)
					((double**)cmd->v2)[ctr][d]=mlist[m]->wrap[d]*space[d]+mlist[m]->pos[d];
				ctr++; }}
		sortVliv((long int*)cmd->v1,cmd->v2,cmd->i1); }

	ctr=0;														// start of code that is run every invokation
	sum=0;
	v1=(long int*)cmd->v1;
	v2=(double**)cmd->v2;
	for(m=0;m<sim->mols->nl[ll];m++)
		if(mlist[m]->ident==i) {
			j=locateVli(v1,mlist[m]->serno,cmd->i1);
			if(j>=0) {
				pos=mlist[m]->pos;
				ctr++;
				if(msddim<0) {
					for(d=0;d<dim;d++) {
						diff=mlist[m]->wrap[d]*space[d]+pos[d]-v2[j][d];
						sum+=diff*diff; }}
				else {
					diff=mlist[m]->wrap[msddim]*space[msddim]+pos[msddim]-v2[j][msddim];
					sum+=diff*diff; }}}
	fprintf(fptr,"%lf %lf\n",sim->time,sum/ctr);

	return 0; }



/****************** internal routines *******************/

int insideecoli(double *pos,double *ofst,double rad,double length) {
	double dist;

	dist=(pos[1]-ofst[1])*(pos[1]-ofst[1])+(pos[2]-ofst[2])*(pos[2]-ofst[2]);
	if(pos[0]-ofst[0]<rad) dist+=(pos[0]-ofst[0]-rad)*(pos[0]-ofst[0]-rad);
	else if(pos[0]-ofst[0]>length-rad) dist+=(pos[0]-ofst[0]-length+rad)*(pos[0]-ofst[0]-length+rad);
	return dist<rad*rad; }


void putinecoli(double *pos,double *ofst,double rad,double length) {
	double dist;

	dist=(pos[1]-ofst[1])*(pos[1]-ofst[1])+(pos[2]-ofst[2])*(pos[2]-ofst[2]);
	if(pos[0]-ofst[0]<rad) {
		dist+=(pos[0]-ofst[0]-rad)*(pos[0]-ofst[0]-rad);
		dist=sqrt(rad*rad/dist);
		pos[0]=ofst[0]+rad+dist*(pos[0]-ofst[0]-rad); }
	else if(pos[0]-ofst[0]>length-rad) {
		dist+=(pos[0]-ofst[0]-length+rad)*(pos[0]-ofst[0]-length+rad);
		dist=sqrt(rad*rad/dist);
		pos[0]=ofst[0]+length-rad+dist*(pos[0]-ofst[0]-length+rad); }
	else
		dist=sqrt(rad*rad/dist);
	pos[1]=ofst[1]+dist*(pos[1]-ofst[1]);
	pos[2]=ofst[2]+dist*(pos[2]-ofst[2]);
	return; }


int molinpanels(simptr sim,int ll,int m,int s,char pshape) {
	int p,ps,dim,npnl;
	double *pos;
	panelptr *pnls,pnl;

	if(s<0) {
		for(s=0;s<sim->srfss->nsrf;s++)
			if(molinpanels(sim,ll,m,s,pshape)) return 1;
		return 0; }

	dim=sim->dim;
	if(pshape=='s') ps=2;
	else return 0;
	pnls=sim->srfss->srflist[s]->panels[ps];
	npnl=sim->srfss->srflist[s]->npanel[ps];
	pos=sim->mols->live[ll][m]->pos;
	if(pshape=='s') {
		for(p=0;p<npnl;p++) {
			pnl=pnls[p];
			if(Geo_PtInSphere(pos,pnl->point[0],pnl->point[1][0],dim)) return 1; }}
	return 0; }