// Simulation and analysis of Spatiotemporal Component of Cl- Gradients in a isolated dendrite

// Define spatial properties of dendrite
  ANZAHL_NODES = 103
  MAX_SYNAPSES = 100

// Determination Parameters GABA
  ANZAHL_GABA = 1        //Anzahl der synaptischen Pulse
  G_GABA = 0.000789     //This Value was pultiplid for the analysis 
  P_GABA = 0.0
  DECAY_GABA = 37  
  ONSET_PULSE = 10  

//----- Define run parameters ---------------------
  tstop = 4000   // Duration 4000
  v_init = -60   // Initial voltage
  dt = 0.5      // Step Interval in ms
  lenghtoutputvec = 6000           //6000     // Number of Lines for output (< 32000 for Excel-Figures)



//----- Insert synapses -------------------------
// Determination of Synapses
  objref gabasyn[MAX_SYNAPSES]
  for i=0, MAX_SYNAPSES-1 {  
  dend {
      printf("%g, ", (i+1)/(ANZAHL_NODES))
      gabasyn[i] = new gaba((i+1)/(ANZAHL_NODES))
      gabasyn[i].tau1 = 0.1 
      gabasyn[i].tau2 = DECAY_GABA
      gabasyn[i].P = P_GABA
    }
  }

// Definition of synaptic Stimuli
  objref stimGABApuls[MAX_SYNAPSES] //Pulssequenz GABA
  for i=0, MAX_SYNAPSES-1 {   
     stimGABApuls[i] = new NetStim((i+1)/(ANZAHL_NODES))
     stimGABApuls[i].number = ANZAHL_GABA
     stimGABApuls[i].start = ONSET_PULSE
  }

// Linkage of synaptic Inputs
  objref synpulsegaba[MAX_SYNAPSES]



// ------------Procedures and Functions -------------------------------
// --------------------------------------------------------------------

// Function MakeShort ---------------------------------------//
// Inputs: $1 Objref to Inputvector                          //
//         $2 Objref to Outoutvector                         //
//         lenoutvec  desired lendth of Outputvector         //
//                                                           //
// Reduce Inputvec to Outputvev by averaging n elements      //
// n (reducing factor) = floor(Inputvec.size() / lenoutvec)  //
// ----------------------------------------------------------//

obfunc MakeShort() {local i, n

  n = int($o1.size()/$3)
  $o2.resize($3)
  for i=0, $3-1 {
    $o2.x[i] = $o1.mean(i*n, (i+1)*n-1)
  }
  return $o2
}   //  End of function


// ---------Definition of Output Vectors and File Output --------------
// --------------------------------------------------------------------

//-- Define  ------
  objref timevec, shorttimevec
  objref voltvec[ANZAHL_NODES], shortvoltvec, voltvec_aver
  objref clivec[ANZAHL_NODES], shortclivec, clivec_aver
  objref Outmatrix  //0=time, SYNAPSES=Cl-, MAX_SYNAPSES+SYNAPSES = v

  strdef OutFileName                       // Name of File Output
  objref OutFile 

//--- Assign -------
  timevec = new Vector()
  shorttimevec = new Vector()
  for i = 0, ANZAHL_NODES-1 {    
     voltvec[i] = new Vector()
     clivec[i] = new Vector()
  }
  clivec_aver = new Vector()
  voltvec_aver = new Vector()
  shortclivec = new Vector()
  shortvoltvec = new Vector()
  Outmatrix = new Matrix()


//-- Link Output Vectors ----------------
  timevec.record(&t)                          // Time vector
  for i = 0, ANZAHL_NODES-1 {                 // Generate Vektor for each node
    voltvec[i].record(&dend.v(i/ANZAHL_NODES))    
    clivec[i].record(&dend.cli(i/ANZAHL_NODES))
  }

// --- Simulation starts here

SYNAPSES = 1
printf(" --------Loop Starts with  %g Synapses of %g ------- \n", SYNAPSES, MAX_SYNAPSES)
while (SYNAPSES <= MAX_SYNAPSES) {

  
  for i = 0, MAX_SYNAPSES-1 {                          // remove conductivity from all synapses 
    dend {
        synpulsegaba[i] = new NetCon(stimGABApuls[i], gabasyn[i], 0, 0, 0)
    }  
  }

  printf(" ---- %g Synapses of %g ------- \n", SYNAPSES, MAX_SYNAPSES)

  for i = 1, SYNAPSES {                                          // activate the right synaspes
    synapse_loc =  int(i/(SYNAPSES+1)*MAX_SYNAPSES)       //in percent for localization
    printf(" - %g - %g - %g  \n", SYNAPSES, i, synapse_loc)
    // Assign Value to Synapse
    dend {
        synpulsegaba[synapse_loc-1] = new NetCon(stimGABApuls[synapse_loc-1], gabasyn[synapse_loc-1], 0, 0, G_GABA)
    }
  }


  // Run Simulation --------------------------------------------------------
     run()
  
  // Put Data in Output Matrix --------------------------------------------
  // ---- => shrink the parameters to output-size before ---------------------
   MakeShort(timevec, shorttimevec, lenghtoutputvec)   
   Outmatrix.resize(shorttimevec.size(), 1 + 2 * MAX_SYNAPSES)
   Outmatrix.setcol(0, shorttimevec)
   
   // Calculate average [Cli] over all Nodes -------
   clivec_aver.resize(clivec[0].size())
   voltvec_aver.resize(voltvec[0].size())
   clivec_aver.mul(0)                      // empty vector
   voltvec_aver.mul(0) 
   for i=0, MAX_SYNAPSES-1 {
        clivec_aver.add(clivec[i])
        voltvec_aver.add(voltvec[i])
    } 
    clivec_aver.div(MAX_SYNAPSES)
    voltvec_aver.div(MAX_SYNAPSES)

   printf("%g Synapses of %g => [Cl-]i = %g mM, Em = %g mV \n", SYNAPSES, MAX_SYNAPSES, clivec_aver.min(), voltvec_aver.max())

   // Put Vectors in Outmatrix 
   MakeShort(clivec_aver, shortclivec, lenghtoutputvec)   
   MakeShort(voltvec_aver, shortvoltvec, lenghtoutputvec)
   Outmatrix.setcol(SYNAPSES, shortclivec)
   Outmatrix.setcol(MAX_SYNAPSES + SYNAPSES, shortvoltvec)
  
   // Goto next SYNAPSES
     SYNAPSES+=1
}  // end of SYNAPSES _ LOOP


// Save the Data --------------------------------------------------------------------
  OutFile = new File()
  sprint(OutFileName, "Result_Isolated_Dendrite_Var_nGABA.asc")
  OutFile.wopen(OutFileName) 
  Outmatrix.fprint(OutFile, "\t%g")     
  OutFile.close