load_file("nrngui.hoc")
cvode_active(0) //1-time step variabile 0-time step fisso
opt= 0	    //Stampa files di corrente per ogni frames
                //0=no,1=yes

tstop=30	    // Durata della simulazione
dist=1

xopen("geoc8076.hoc")       //Lettura file morfologia      
numaxon=1	    // axon
numsoma=1	    // soma
numbasal=77		// dendrites
numapical=48	// apical_dendrite=obliqui
                // Costituiscono la max parte della superficie apicale e
				// sono sede principale dei bersagli dei assoni collaterali 
				// di Schaff delle CA3.
numtrunk=51 	// user 5

//***************************************************************************** 
// PROPRIETA' ELETTRICHE (Non modificare!!)
//***************************************************************************** 
Rm =    28000.		  // Ohm*cm**2
RmDend = Rm/1.
RmSoma = Rm
RmAx = Rm											

Cm    = 1.		  // uF/cm**2
CmSoma= Cm
CmAx  = Cm
CmDend = Cm*1.

RaAll= 150.		  // Ohm*cm Resistenza intracellulare
RaSoma=150.  	  // Ohm*cm
RaAx = 50.	      // Ohm*cm Resistenza nell' assone

Vrest = -65.	  // mV
dt = 0.013
gna =  .025
AXONM = 5.
gkdr = 0.01
celsius = 35.0    // Gradi Celsius
KMULT =  0.03
KMULTP = 0.03
ghd=0.00005
nash=0.
//*****************************************************************************

//**************************************************************************
// DETERMINAZIONE NUMERO DI COMPARTIMENTI
//**************************************************************************
// In questo caso il numero di compartimenti è uguale al numero dei segmenti
// col quale è stato discretizzato il neurone e ricostruita la morfologia.
// E' stato visto che la lunghezza tipica di cavo per le proprietà biofisiche
// del presente neurone è molto più grande delle dimensioni tipiche di ciascun
// compartimento, così la differenza di potenziale tra adiacenti compartimenti
// è molto piccola. Ciò ci interessa perchè vogliamo calcolare il campo
// magnetico del neurone con un grado di precisione molto alto, ed abbiamo
// bisogno di un numero molto elevato di segmentazioni del neurone.

  forall {nseg=n3d()-1} // Assegna un numero di compartimenti entro
                        // ciascuna sezione pari al numero di segmenti
					    // in cui è stato sezionato.

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

//************************************************************************ 
//INSERIMENTO COMPONENTI PASSIVE
forsec "axon" {insert pas e_pas=Vrest g_pas = 1./RmAx Ra=RaAx cm=CmAx}
forsec "soma" {insert pas e_pas=Vrest g_pas = 1./RmSoma Ra=RaSoma cm=CmSoma}
forsec "dendrite"{insert pas e_pas=Vrest g_pas = 1./RmDend Ra=RaAll cm=CmDend}
forsec "user5" {insert pas e_pas=Vrest g_pas = 1./RmDend Ra=RaAll cm=CmDend}
//************************************************************************

objref b,c, stim, outfile,outfile1
objref sref, blist[numtrunk]
strdef filename, dend, trunk ,fil,filename1,fil1

access soma


distance()  // Fissa il riferimento a x=0.5 nella sezione corrente
            // in questo caso la sezione corrente è il soma. La
			// distanza dal soma è importante perchè a seconda di
			// essa varia la densità dei vari canali ionici.

xopen("synapses.hoc")  //Connessioni sinaptiche

//***************************************************************************** 
// INSERIMENTO COMPONENTI ATTIVE
//***************************************************************************** 

forsec "axon" {   
                insert nax gbar_nax=gna * AXONM	sh_nax=nash
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP*1
 		insert ds
}

forsec "soma" {   
		insert hd ghdbar_hd=ghd	vhalfl_hd=-73
                insert na3 ar_na3=1 sh_na3=nash gbar_na3=gna
                insert kdr gkdrbar_kdr=gkdr
                insert kap gkabar_kap = KMULTP
 		insert ds
}

for i=0, numbasal-1 dendrite[i] {
 		insert ds
		if (diam>0.35) {factor=1} else {factor=1}
		insert hd ghdbar_hd=ghd vhalfl_hd=-73
                insert na3 ar_na3=1 gbar_na3=gna*factor sh_na3=nash
                insert kdr gkdrbar_kdr=gkdr*factor
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0

		for (x) if (x>0 && x<1) { xdist = distance(x)
                	ghdbar_hd(x) = factor*ghd*(1+3*xdist/100)
                		if (xdist > 100){
					vhalfl_hd=-81
                        		gkabar_kad(x) = factor*KMULT*(1+xdist/100)
                			} else {
					vhalfl_hd=-73
                        		gkabar_kap(x) = factor*KMULTP*(1+xdist/100)
               				}
		}
}
                
forsec "apical_dendrite" {
	insert ds
	if (diam>0.35) {factor=1} else {factor=1}
		insert hd ghdbar_hd=ghd
                insert na3 ar_na3=1 gbar_na3=gna*factor sh_na3=nash
                insert kdr gkdrbar_kdr=gkdr*factor
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0

		for (x) if (x>0 && x<1) { xdist = distance(x)
                	ghdbar_hd(x) = factor*ghd*(1+3*xdist/100)
                		if (xdist > 100){
					vhalfl_hd=-81
                        		gkabar_kad(x) = factor*KMULT*(1+xdist/100)
                			} else {
					vhalfl_hd=-73
                        		gkabar_kap(x) = factor*KMULTP*(1+xdist/100)
               				}
		}
}

forsec "user5" {
	insert ds
		insert hd ghdbar_hd=ghd
                insert na3 ar_na3=1 gbar_na3=gna sh_na3=nash
                insert kdr gkdrbar_kdr=gkdr
		insert kap gkabar_kap=0
		insert kad gkabar_kad=0

		for (x) if (x>0 && x<1) { xdist = distance(x)
                	ghdbar_hd(x) = ghd*(1+3*xdist/100)
                		if (xdist > 100){
					vhalfl_hd=-81
                        		gkabar_kad(x) = KMULT*(1+xdist/100)
                			} else {
					vhalfl_hd=-73
                        		gkabar_kap(x) = KMULTP*(1+xdist/100)
               				}
		}
}
//************************************************************** 

//*****************************************************
// GRAFICA //
//******** Per creare una box *************************
objref vbox,vbox1
objref g,h,n,p,ll

vbox = new VBox() 
vbox.intercept(1) 
g = new Graph()
g.size(0,tstop,-70,30)

g.addvar("apical_dendrite[12].v(rel)",1,1,0.5,1,2)
g.addvar("apical_dendrite[14].v(rel)",2,1,0.5,0.98,2) 
g.addvar("apical_dendrite[30].v(rel)",3,1,0.5,0.96,2) 

g.xaxis(1)

h = new Graph()
h.size(0,tstop,-70,30)

h.addvar("soma.v(0.5)",1,1,0.5,1,2)

h.xaxis(1)

xpanel("",0)
xbutton("Start-Restart","runm()")
xpanel()  
vbox.intercept(0) 
vbox.map("Potenziali (mV)",427,98,322,540)	   //b = new VBox()

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

vbox1 = new VBox() 
vbox1.intercept(1) 

n = new Graph()	   // Dipoli
n.size(0,tstop,-3.,3.)
n.addexpr("dipx",1,1,0.5,1,2)
n.addexpr("dipy",2,1,0.5,0.98,2)
n.addexpr("dipz",3,1,0.5,0.96,2)
n.addexpr("diptot",4,1,0.5,0.94,2)
n.xaxis(1)

vbox1.intercept(0) 
vbox1.map("Fig.8a ECD (pA*m)",-2,98,322,510.3)

p = new PlotShape()
p.exec_menu("Shape Plot")
p.size(-245.834,295.834,-671.282,341.282) 
p.label(0,0,"geoc8076",2,1,0,0,1)
p.variable("v")
p.show(1)


//***************************************************************************** 
// INIZIALIZZAZIONE SIMULAZIONE
//***************************************************************************** 
			
proc init() {
	t=0
        conto=0
        forall {
        v=Vrest
        if (ismembrane("nax") || ismembrane("na3")) {ena=55}
        if (ismembrane("kdr") || ismembrane("kap") || ismembrane("kad")) {ek=-90}
        if (ismembrane("hd") ) {ehd_hd=-30}
	}
	finitialize(Vrest)
      fcurrent()

      forall {
	for (x) {
	if (ismembrane("na3")||ismembrane("nax")){e_pas(x)=v(x)+(ina(x)+ik(x))/g_pas(x)}
	if (ismembrane("hd")) {e_pas(x)=e_pas(x)+i_hd(x)/g_pas(x)}
		}
	}
	cvode.re_init()
	cvode.event(tstop)
	access soma
	g.begin()
      h.begin()
      n.begin()

}

    proc advance() {
    fadvance()

//***********************************************************************
// CALCOLO DEL DIPOLO DI CORRENTE (Q)
//***********************************************************************
// Da Murakami, Hirose, Okada, (J. Physiol. (2003) 553,3,pp. 975-985)
// "Contr. of ionic curr. to MEG and EEG signals generated by guinea 
//  CA3 slices"
//***********************************************************************

      dipx=0.
	dipolox=0.
	dipy=0.
	dipoloy=0.
	dipz=0.
	dipoloz=0.
	diptot=0.
	
	forall for i=0, n3d()-2 {
	
	resist=ri((i+0.)/(n3d()-1))	  //  MOhm
	difftens=( v((i+0.)/(n3d()-1))-v((i+1)/(n3d()-1)) )	 //Diff. Tens. Adiac. Comp.
	curr=difftens/resist          // nA
	
	dipolox=(x3d(i)-x3d(i+1))*curr*0.001	 // in pA*m
	dipoloy=(y3d(i)-y3d(i+1))*curr*0.001
	dipoloz=(z3d(i)-z3d(i+1))*curr*0.001	
	dipx=dipx+dipolox
	dipz=dipz+dipoloz
	dipy=dipy+dipoloy						

	diptot=sqrt(dipx*dipx+dipy*dipy+dipz*dipz) // in pA*m

}
	g.plot(t)
	h.plot(t)
	n.plot(t)
      g.flush()
	h.flush()
	n.flush()
	p.flush()
      doNotify()
}

proc runm() { 
run()
}