// processes and functions 

//fixing nseg number of compartments for the spatial discretization

// these are reasonable values for most models
freq = 100      // Hz, frequency at which AC length constant will be computed
d_lambda = 0.1

func lambda_f() { local i, x1, x2, d1, d2, lam
        if (n3d() < 2) {
                return 1e5*sqrt(diam/(4*PI*$1*Ra*cm))
        }
// above was too inaccurate with large variation in 3d diameter
// so now we use all 3-d points to get a better approximate lambda
        x1 = arc3d(0)
        d1 = diam3d(0)
        lam = 0
        for i=1, n3d()-1 {
                x2 = arc3d(i)
                d2 = diam3d(i)
                lam += (x2 - x1)/sqrt(d1 + d2)
                x1 = x2   d1 = d2
        }
        //  length of the section in units of lambda
        lam *= sqrt(2) * 1e-5*sqrt(4*PI*$1*Ra*cm)

        return L/lam
}

//initializing the simulation

proc init(){

	forall{
		for(x,0){
			Ra = user_Ra_b
			v = 0
			//cm = user_cm_b
			//g_pas = 1/user_Rm_b
			e_pas = user_e_pas
			v_init = 0
		}
	}
	
	forsec a_source_blue[0].all{
		for(x,0){
			g_pas = 1/user_Rm_b
			cm = user_cm_b
		}
	}
	
	forsec a_source_red[0].all{
		for(x,0){
			g_pas = 1/user_Rm_r
			cm = user_cm_r
		}
	}
	
	ntarget = 10
	
	for(i=0; i<ntarget; i+=1){
		forsec a_target_blue[i].all{
			for(x,0){
				g_pas = 1/user_Rm_b
				cm = user_cm_b
			}
		}
		forsec a_target_red[i].all{
			for(x,0){
				g_pas = 1/user_Rm_r
				cm = user_cm_r
			}
		}
	}
	
	for(i=0; i<ntarget*2; i+=1){
	
		elecsyn_NetConn_source_blue_target_blue_gapCond_A[i].weight = gapWeight 
		elecsyn_NetConn_source_blue_target_blue_gapCond_B[i].weight = gapWeight 
		
		elecsyn_NetConn_source_red_target_red_gapCond_A[i].weight = gapWeight
		elecsyn_NetConn_source_red_target_red_gapCond_B[i].weight = gapWeight		
	}
	
	for(i=0; i<4; i+=1){
		elecsyn_NetConn_source_blue_source_red_gapCond_A[i].weight = gapWeight
		elecsyn_NetConn_source_blue_source_red_gapCond_B[i].weight = gapWeight
	}
	
	finitialize(v_init)
    if (cvode.active()) {
      cvode.re_init()
    } else {
      fcurrent()
    }
    frecord_init()
}

//Setting dt

proc setdt() {local Dt, dtnew
	if (using_cvode_) return
	Dt = 1/steps_per_ms
	nstep_steprun = int(Dt/dt)
	if (nstep_steprun == 0) {
		nstep_steprun = 1
	}
	dtnew = Dt/nstep_steprun
	if (abs(dt*nstep_steprun*steps_per_ms - 1) > 1e-6) {
		print "Changed dt"
		dt = dtnew
	}
}

//Setting the screen update interval

proc continuerun() {local rt, rtstart, ts
	realtime = 0  rt = screen_update_invl  rtstart = startsw()
	eventcount=0
	eventslow=1
	stoprun = 0
	if (using_cvode_) {
		cvode.event($1)
		ts = $1
		if (cvode.use_local_dt) {
			cvode.solve(ts)
			flushPlot()
			realtime = startsw() - rtstart
			return
		}
	}else{
		ts = $1 - dt/2
	}
	while(t < ts && stoprun == 0) {
		step()
		realtime = startsw() - rtstart
		if (realtime >= rt) {
//			if (!stdrun_quiet) fastflushPlot()
			screen_update()
			//really compute for at least screen_update_invl
			realtime = startsw() - rtstart
			rt = realtime + screen_update_invl
		}
	}
	if (using_cvode_ && stoprun == 0) { // handle the "tstop" event
		step() // so all recordings take place at tstop
	}
	flushPlot()
	realtime = startsw() - rtstart
}


//Possibility to choose between the integration methods

proc cvode_act(){
	if (isCVodeAct){
		cvode.active(1)
		
		print "Using variable time step integration method - faster simulation"
	} else{
		cvode.active(0)
		
		print "Using fix time step integration method - slower simulation"
	}
}

proc injSoma_b(){

	isStim1 = 1
	isStim2 = 0

	stim1.amp = 0.2
	stim2.amp = 0.05
	stim3.amp = 0
	stim4.amp = 0
}

proc injSoma_r(){

	isStim1 = 0
	isStim2 = 1
/*
	access soma_r{
		stim3.loc(0.5)					
		stim4.loc(0.5)				
	}
*/
	stim1.amp = 0
	stim2.amp = 0
	stim3.amp = 0.2
	stim4.amp = 0.05
}