mesh.hoc

// Mesh for NEURON
// (c) Maciej Lazarewicz (mlazarew@gmu.edu) (2001) v1.1

strdef maxsec,maxdifsec
objref imp1

// This routine mash the cell that each segment has length smaller then $2*(length constant lambda_f) 
// for frequency $1 (Hines M.L., Carnavale N.T. Neuroscientist, 7(2):123-135, 2001 )

func lambdaf() {
	if($2>0) {	
		return 1e5*sqrt(diam($1)/(4*PI*$2*Ra*cm))
	} else {
		return 1e2*sqrt(diam($1)/(4*Ra*g_pas))
	}
}

//If you use pt3d structures, before you want to calculate lambda_100, you have to set up nseg for desire value, becauce diam can depend on nseg
func mesh() {
	nr=0
	max=0
	forall nseg=1
	finitialize()
	forall {
		while(1) {
		  maxl=0
		  finitialize()
		  for (x) {
                  	el=L/(nseg*lambdaf(x,$1)) 
                  	if(maxl<el) {
				maxl=el
			}  
                  }
		  if(maxl<=$2)break
		  nseg+=1
		}
  		if(maxl>max) {
			max=maxl
			maxsec=secname()
		}
		nr+=nseg
	}
	finitialize()
	print "Max l100 =",maxsec,max,"#comp:",nr
	return nr
}

func ratio_mesh() {
	imp1=new Impedance()
	nr=0
	maxdif=0
	forall nseg=1
	finitialize()
	forall {
		while(1) {
		  if(issection("soma"))break
		  maxl=0
		  finitialize()
		  imp1.loc(0)
		  imp1.compute($1)
		  tra=imp1.ratio(1/nseg)
		  if(abs(1-tra)<$2)break
		  nseg+=1
		}
  		if(maxdif<tra) {
			maxdif=tra
			maxdifsec=secname()
		}
		nr+=nseg
		print secname(), " ",nseg
	}
	print "Maxdif =",maxdifsec,maxdif,"#comp:",nr
	finitialize()
	return nr
}

func transfer_mesh_2() {
	imp1=new Impedance()
	nr=0
	maxdif=0
	forall nseg=1
	finitialize()
	forall {
		prev=0
		while(1) {
		  xopen("Hoc/utilities.hoc")maxl=0
		  finitialize()
		  imp1.loc(0.5)
		  imp1.compute($1)
		  tra=imp1.input(0.5)
		  dif=abs(tra-prev)
//		  print nseg,tra
		  if(dif<$2 && prev!=0)break
		  prev=tra
		  nseg+=2
		}
  		if(maxdif<dif) {
			maxdif=dif
			maxdifsec=secname()
		}
		nr+=nseg
		print secname(), " ",nseg
	}
	print "Maxdif =",maxdifsec,maxdif,"#comp:",nr
	finitialize()
	return nr
}
objref diff,maxx

proc make_table() {
	diff=new Vector(10000)
	maxx=new Vector(10000)
	forall nseg=9
	finitialize()
	init()
	print "Start:",t,tstop
	while(t<tstop){
		step()
		ii=0
		forall {
			diff.x[ii]=abs(v(0)-v(1))
			if(diff.x[ii]>maxx.x[ii])maxx.x[ii]=diff.x[ii]
			ii+=1
		}
	}
	print "Stop"
	ii=0
	forall {
		print maxx.x[ii]
		ii+=1
	}
}

func drop_mesh() {
	nr=0
	ii=0
	forall {
		nseg=int(maxx.x[ii]/$1)+1
		nr+=nseg
		ii+=1
	}
	print "# comp: ",nr
	
	finitialize()
	return nr-soma.nseg+1
}


func length_mesh() { local nr
	print "I started meshing"
	nr=0
	forall {
		nseg=int(L/$2)+1                  	
		nr+=nseg
	}
	finitialize()
	print "#comp:",nr
	return nr
}
objref onseg
func fast_mesh() { local nr,max,el,i
	print "I started meshing"
	nr=0
	onseg=new Vector(10000)
	i=0
	max=0
	forall { 
		onseg.x[i]=nseg 
		i+=1
	}
	forall nseg=1
	finitialize()
	i=0
	forall {
		el=L/(lambdaf(0.5,$1))
		nseg=int(el/$2)+1                  	
		if(max<el/nseg) {
			max=el/nseg
			maxsec=secname()
		}
//		if(onseg.x[i]!=nseg)printf("I changed segmend %s from %d to %d\n",secname(),onseg.x[i],nseg)
		i+=1
		nr+=nseg
	}
	finitialize()
	print "Max l100 =",maxsec,max,"#comp:",nr
	return nr
}

func fast_mesh_one() {
	nr=0
	max=0
	nseg=1
	finitialize()
		el=L/(lambdaf(0.5,$1))
		nseg=int(el/$2)+1                  	
		if(max<el/nseg) {
			max=el/nseg
			maxsec=secname()
		}
		nr+=nseg	
	finitialize()
	print "Max l100 =",maxsec,max,"#comp:",nr
	return nr
}

func uniform_mesh() {

	forall nseg=$1
	soma.nseg=1
	finitialize()		
	return 210*$1+1
}

objref inj2,imp,eir,elc,f22
strdef n1

proc comp() {local th11,i
	
	eir=new Vector(1000)
	elc=new Vector(1000)

	f22=new File()

	sprint(n1,"%s.dat",$s1)
	f22.wopen(n1)

	xl=0.5
  inj2 = new IClamp(0.5)
  inj2.del=10  //10
  inj2.dur=3  //3
  inj2.amp=1
	soma {		
		imp=new Impedance()

		max_nr_comp=6000
		i=1
		make_table()
		inj2.amp=0
		while(elc.x[i-1]<max_nr_comp) {
			th11=100/i
			print "Th: ",th11
//			nr=0
//			forall {
//				nseg=i
//				nr+=nseg
//			}
//			elc.x[i]=nr+1-i
//			print elc.x[i]
			elc.x[i]=drop_mesh(th11)	
//			elc.x[i]=fast_mesh(0,th11)	
//			forall nseg=3^(i-1)
//		        elc.x[i]=0
//			forall elc.x[i]+=nseg
			soma.nseg=1
			soma distance()
//			ins_active()

//			elc.x[i]=uniform_mesh(i)
//			elc.x[i]=fast_mesh_one(100,th11)	
//			elc.x[i]=ratio_mesh(100,th11)	
//			elc.x[i]=transfer_mesh_2(100,th11)	

			init()
			f22.printf("%d %d ",i,elc.x[i])
			for(k=0;k<=20;k+=1) {
			soma imp.loc(xl)
			soma	imp.compute(10*k)
			soma	eir.x[i]=imp.input(xl)

				f22.printf("%10.1000g ",eir.x[i])
			}
			f22.printf("\n")
			i+=1
		}
		f22.close()
	}
}

objref w11,w22,file123
strdef name
proc comp2() {local th11,i,j
	file123=new File()
	for(i=0;i<1;i+=1) {
		if(i>0)th11=1/i else th11=0.04
		print "Th: ",th11
		el1=fast_mesh(100,th11)
		soma.nseg=1
//		el1=uniform_mesh(i)
//		el1=fast_mesh_one(100,th11)	
//		el1=ratio_mesh(100,th11)	
//		el1=transfer_mesh_2(100,th11)	

 		w11=new Vector(30000)
		w22=new Vector(30000)

		cvode.record(&soma.v(0.5),w22,w11)
		soma distance()
		ins_active()
		init()
		while(t<tstop) step()
		sprint(name,"Data1/%s_%0.6d_%0.5d.dat",$s1,el1,i)
		file123.wopen(name)
		for(j=0;j<w11.size();j+=1)file123.printf("%10.100g %10.100g\n",w11.x[j],w22.x[j])
		file123.close()
		
	}
}

objectvar rvp,yrvp,xrvp,file12,fdx,fdy

proc save_spatial() {
   init()

   while(t<tstop) {
	   step()
	   file12 = new File()
	   xrvp   = new Vector()
     	   yrvp   = new Vector()
	   fdx    = new Vector()
     	   fdy 	  = new Vector()
	   fdx    = fdx.indgen(-265,804,1)
	
           rvp = new RangeVarPlot("v")
dend8[10]  rvp.begin(1)
dend3[83]  rvp.end(1)
           rvp.to_vector(yrvp,xrvp)
	   fdy = yrvp.c.interpolate(fdx,xrvp)
           sprint(name,"Images/%f.ms",t)
           file12.wopen(name)
   	   for(j=0;j<fdy.size();j+=1)file12.printf("%g %10.100g %10.100g\n",t,fdx.x[j],fdy.x[j])
	   file12.close()
  }
}

proc disc() {
forall {
	if ((distance(0)>100 || distance(1)>100 ) && !issection("axon")) for (x) { 
disconnect()
//		v=Vrest
//		ina=0
//		ik=0
//		ica=0 
	}
}
}