load_file("nrngui.hoc")
create soma[1]
objectvar st
objectvar sh
strdef str1

tstop = 20
steps_per_ms = 40
ra     = 200			
rm     = 60000
c_m       = 0.75

v_init    = -65
celsius   = 30


proc setpassive() {
	forall {
	  	insert pas
	  	cm = c_m 
	  	e_pas = v_init
		Ra=ra
		g_pas=1/rm
	}
}

proc load_3dcell() {
	xopen(str1)
	setpassive()
	forall { 
		soma area(0.5)
		nseg = (int((L/(0.1*lambda_f(100))+.9)/2)*2 + 1)
		totcomp=totcomp+nseg
	}
	setpassive()
	access soma[0]
}	

/*****************************************************************/

objref netlist, s, ampasyn, f1, DEND, sapamp, somavec, sampvec
strdef  str2

proc find_epsp_amplitudes() {
	nmdaamparat=0.25
	atrise=2.0
	atau=3
	ntrise=5
	ntau=50

    netlist=new List()

	
    ropen("neuron.que")
	f1=new File()
	nneu=fscan()
	tstop=40
    Gstart=5e-5
    Gend=1
	while (nneu > 0){
    	getstr(str1,1)
		load_3dcell(str1)
		sprint(str2, "%s.nor",str1)
		wopen(str2)
		
		soma[2] s = new NetStim(0.5)
    	s.interval=1   // ms (mean) time between spikes
    	s.number=1     // (average) number of spikes
    	s.start=2   // ms (mean) start time of first spike
    	s.noise=0      // range 0 to 1. Fractional randomness.
		
		DEND = new SectionList()
		forsec "basal" { 
			DEND.append()
		}
		forsec "apical"{
			DEND.append()
		}
		access soma[2]
		distance()

		forsec DEND {
			for (x) {
    			G=Gstart
				ampasyn=new AMPA(x)
				ampasyn.TRise=atrise
				ampasyn.tau=atau
    			netlist.append(new NetCon(s,ampasyn,1,0,0))

    			while(G<=Gend){
					ampasyn.gmax=G
        			update_init()
            		
					MAX=-70
					DENDVMAX=-70
					while (t <tstop){
                    	fadvance()
						if(soma[2].v(0.5)>MAX) {
							MAX=soma[2].v(0.5)
						}	
						if (v(x)>DENDVMAX) {
							DENDVMAX=v(x)
						}	
        			}
	
					if(MAX>-64.8) { // somatic EPSP amplitude crossed 0.2 mV
						print secname(), " ", x, " ", distance(x), " ", 65+DENDVMAX, " ", G, " ",  65+MAX
            			fprint("%f\t%f\t%f\n", distance(x), 65+DENDVMAX, G)
        				update_init()
						break
					}	
        			G=G+1e-5
    			}
				if(G>=Gend){ // Even after reaching Gend, Soma EPSP amplitude is not 0.2 mV
						print secname(), " ", x, " ", distance(x), " ", 65+DENDVMAX, " ", "-1", " ", 65+MAX
            		fprint("%f\t%f\t-1\n",distance(x),65+DENDVMAX)
				}	
					
			}
		}	
    	wopen()
		nneu=nneu-1
	}
}

/*****************************************************************/

proc update_init(){
    finitialize(v_init)
    fcurrent()
}	
 
/*****************************************************************/

find_epsp_amplitudes()