/*---------------------------------------------
                   MERGING PROCEDURES
  ---------------------------------------------*/

// Useful variables
objectvar clusterVectApp



// merging procedure for soma sections
proc mergingSoma(){local nsect,sec,sec0,totsurf,totri,newL,newRa,newri,Ra1,Ra2,L1,L2,diam1,diam2,ri1,ri2,newdiam1
     clusterVectApp = $o1
     nsect = clusterVectApp.size()-1
     sec0 = clusterVectApp.x(0)
     totsurf = 0
     totri = 0
     newL = 0
     newRa = 0
     newri = 0
     for ii=0,nsect {
         sec = clusterVectApp.x(ii)
         totsurf = totsurf + secSurf.x(sec)
         totri = totri + secri.x(sec)
         newL = newL + secL.x(sec)
         newRa = newRa + secRa.x(sec)
         newri = newri + secri.x(sec)
         secpathri.set(sec,totri)
     }
     newRa = newRa/(nsect+1)
     if (PRESERVINGsomaMETHOD==0){
        newdiam = sqrt(newRa*newL*4/newri/PI/100)
     }else if (PRESERVINGsomaMETHOD==1){
        newdiam = totsurf/(newL*PI)
     }
     
     eqL.set(0,newL)
     eqdiam.set(0,newdiam)
     eqri2.set(0,newri)
     orSurf.set(0,totsurf)
     orSurf1.set(0,totsurf)
     eqSurf.set(0,newL*PI*newdiam)
     surfFact.set(0,surfFactVect.x(0)*orSurf.x(0)/eqSurf.x(0))
     surfFactRmCm.set(0,surfFactVect.x(0)*orSurf.x(0)/eqSurf.x(0))
}

// merging procedure for axon sections
proc mergingAxon(){local nsect,sec,sec0,totsurf,totri,newL,newRa,newri,Ra1,Ra2,L1,L2,diam1,diam2,ri1,ri2,newdiam1
     clusterVectApp = $o1
     nsect = clusterVectApp.size()-1
     sec0 = clusterVectApp.x(0)
     totsurf = 0
     totri = 0
     newL = 0
     newRa = 0
     newri = 0
     for ii=0,nsect {
         sec = clusterVectApp.x(ii)
         totsurf = totsurf + secSurf.x(sec)
         totri = totri + secri.x(sec)
         newL = newL + secL.x(sec)
         newRa = newRa + secRa.x(sec)
         newri = newri + secri.x(sec)
         secpathri.set(sec,totri)
     }
     newRa = newRa/(nsect+1)
     if (PRESERVINGsomaMETHOD==0){
        newdiam = sqrt(newRa*newL*4/newri/PI/100)
     }else if (PRESERVINGsomaMETHOD==1){
        newdiam = totsurf/(newL*PI)
     }
     
 
     
     eqL.set(1,newL)
     eqdiam.set(1,newdiam)
     eqri2.set(1,newri)
     orSurf.set(1,totsurf)
     orSurf1.set(1,totsurf)
     eqSurf.set(1,newL*PI*newdiam)
     surfFact.set(1,surfFactVect.x(1)*orSurf.x(1)/eqSurf.x(1))
     surfFactRmCm.set(1,surfFactVect.x(1)*orSurf.x(1)/eqSurf.x(1))
}


// merging method for serial sections
//         input:  {Ra1,Ra2,L1,L2,diam1,diam2,ri1,ri2,surf1,surf2,diamSer1,diamSer2,diamPar1,diamPar2}
//         output: {newRa,newL,newdiam,newri,newri2,newdiamSer,newdiamPar}
proc mergingSerialMethod(){local Ra1,Ra2,L1,L2,diam1,diam2,ri1,ri2,surf1,surf2
    {Ra1=$1 Ra2=$2 L1=$3 L2=$4 diam1 =$5 diam2=$6 ri1=$7 ri2=$8 surf1=$9 surf2=$10 diamSer1=$11 diamSer2=$12 diamPar1=$13 diamPar2=$14 }
     
    newRa=(Ra1+Ra2)/2
    if (LENGTHMERGINGserialMETHOD == -1){ //LENGTHMERGINGserialMETHOD = -1 by Destexhe
        newL=L1+L2
    } else if(LENGTHMERGINGserialMETHOD==0) {
		newL=(L1+L2)/2
	} else if(LENGTHMERGINGserialMETHOD==1){
		newL=(L1*diam1+L2*diam2)/(diam1+diam2)
	} else if(LENGTHMERGINGserialMETHOD==2) {
		newL=(L1*surf1+L2*surf2)/(surf1+surf2)	
	}    
    
    newri = ri1+ri2   
    newri2 = newri
    newSurf = max(surf1,surf2)
    newdiam = 1      //init value
    if (PRESERVINGMETHOD==0){             //PRESERVINGMETHOD = 0 DIAMMERGINGserialMETHOD = 0 by Destexhe
       if (DIAMMERGINGserialMETHOD==0){
          newdiam = sqrt(newRa*newL*4/newri/PI/100)
       }else if (DIAMMERGINGserialMETHOD==1 || DIAMMERGINGserialMETHOD==-1){
          newdiamSer = diamSer1+diamSer2
          newdiamPar = diamPar1+diamPar2
       }else if (DIAMMERGINGserialMETHOD==2){
          newdiam = (diam1+diam2)/2
       }
    }else if (PRESERVINGMETHOD==1){
       newdiam = (surf1+surf2)/(newL*PI)
    }    
    
}

// merging method for parallel sections
//         input:  {Ra1,Ra2,L1,L2,diam1,diam2,diamSer1,diamSer2,diamPar1,diamPar2}
//         output: {Ra12,L12,diam12,cross12,surf1,surf2,diam12Ser,diam12Par}
proc mergingParallelMethod(){local Ra1,Ra2,L1,L2,diam1,diam2
    {Ra1=$1 Ra2=$2 L1=$3 L2=$4 diam1 =$5 diam2=$6 diamSer1=$7 diamSer2=$8 diamPar1=$9 diamPar2=$10}
     
	surf1=L1*PI*diam1
	surf2=L2*PI*diam2
    
    if(LENGTHMERGINGparallelMETHOD==-1) {
		L12=L1+L2
	}else if(LENGTHMERGINGparallelMETHOD==0) {
		L12=(L1+L2)/2
	} else if(LENGTHMERGINGparallelMETHOD==1){
		L12=(L1*diam1+L2*diam2)/(diam1+diam2)
	} else if(LENGTHMERGINGparallelMETHOD==2) {
		L12=(L1*surf1+L2*surf2)/(surf1+surf2)	
	} else if(LENGTHMERGINGparallelMETHOD==3) {
		if (L1>L2) {L12=L1} else {L12=L2}
	}
	
    Ra12 = (Ra1+Ra2)/2
    diam12 = 1        //init value
    if (PRESERVINGMETHOD==0){             //PRESERVINGMETHOD = 0 DIAMMERGINGparallelMETHOD=0 by Destexhe

       if (DIAMMERGINGparallelMETHOD==0){
          diam12=sqrt(diam1^2+diam2^2)
       }else if (DIAMMERGINGparallelMETHOD==1 || DIAMMERGINGparallelMETHOD==-1){
          diamSer12 = diamSer1+diamSer2
          diamPar12 = diamPar1+diamPar2
       }       
    }else if (PRESERVINGMETHOD==1){
       diam12 = (surf1+surf2)/(L12*PI)
    }         
    
    cross12=PI*(diam12^2)/4     
}



// merging method for Y_group_sections
//         input:  {Ra1,Ra2,Ra3,L1,L2,L3,diam1,diam2,diam3,ri1,ri2,ri3,diamSer1,diamSer2,diamSer3,diamPar1,diamPar2,diamPar3}
//         output: {newRa,newL,newdiam,newri,newri2}
proc mergingYMethod(){ local Ra1,Ra2,Ra3,L1,L2,L3,diam1,diam2,diam3,ri1,ri2,ri3,ri12,surf3,surf12
    {Ra1=$1 Ra2=$2 Ra3=$3 L1=$4 L2=$5 L3=$6 diam1 =$7 diam2=$8 diam3=$9 ri1=$10 ri2=$11 
     ri3=$12 diamSer1=$13 diamSer2=$14 diamSer3=$15 diamPar1=$16 diamPar2=$17 diamPar3=$18}

    mergingParallelMethod(Ra1,Ra2,L1,L2,diam1,diam2,diamSer1,diamSer2,diamPar1,diamPar2)
	
	ri12 = Ra12*L12/cross12/100
	surf12 = L12*diam12*PI
    surf3 = L3*diam3*PI
    
    
    mergingSerialMethod(Ra12,Ra3,L12,L3,diam12,diam3,ri12,ri3,surf12,surf3,diamSer12,diamSer3,diamPar12,diamPar3)
  
    newri2 = (ri1*ri2/(ri1+ri2))+ri3
    newSurf = max(surf1+surf2,surf3)
}


objectvar weightsVect   // Vector for the i-th cluster, containing the id-sections 

// merging method for all sections after merging parent sections and Y_group_sections
//         input:  {clusterVectApp,index}


proc mergingAllMethod(){local index,nsectOrigin,nsect,sec
    clusterVectApp = $o1
    index = $2
    nsectOrigin = $3
    nsect = clusterVectApp.size()-1
    weightsVect = new Vector(nsect+1,1)
    normalFact = nsect+1
     
    maxL=0
    totsurf=0
    totdiam=0
    eqri2Prod = 1
    eqri2Sum = 0
    surfSum = 0
	for ii=0,nsect {
	    sec = clusterVectApp.x(ii)	
        if(maxL<mrgL.x(sec)) maxL=mrgL.x(sec)
        surf=mrgL.x(sec)*PI*mrgdiam.x(sec)
        //print "mergingAllMethod____:",secNameList.o(sec).s,surf
        totsurf=totsurf+surf
        totdiam=totdiam+mrgdiam.x(sec)
        eqri2Prod = eqri2Prod*mrgri2.x(sec)
        eqri2Sum = eqri2Sum+mrgri2.x(sec)
        surfSum = surfSum+mrgSurf.x(sec)
        
        if(LENGTHMERGINGparallelMETHOD==-1){
              normalFact=nsectOrigin weightsVect.set(ii,1)
        }else if(LENGTHMERGINGparallelMETHOD==0){
              normalFact=nsect+1 weightsVect.set(ii,1)
        }else if(LENGTHMERGINGparallelMETHOD==1){
              normalFact=totdiam weightsVect.set(ii,mrgdiam.x(sec))
        }else if(LENGTHMERGINGparallelMETHOD==2){
              normalFact=totsurf weightsVect.set(ii,surf)
        }else if(LENGTHMERGINGparallelMETHOD==3){
              normalFact=nsect+1 weightsVect.set(ii,1)}        
    }    
   	for ii=0,nsect {
	    sec = clusterVectApp.x(ii)
        if(LENGTHMERGINGparallelMETHOD==3) {eqL.set(index, maxL)
        }else{eqL.set(index, eqL.x(index)+mrgL.x(sec)*weightsVect.x(ii)) }  
         
        eqdiamSer.set(index, eqdiamSer.x(index)+ mrgdiamSer.x(sec))    
        eqdiamPar.set(index, eqdiamPar.x(index)+ mrgdiamPar.x(sec))    
        
        if(RESCALEMOD==0){                  
             eqdiam.set(index, eqdiam.x(index)+mrgdiam.x(sec)^2)
        }else if (RESCALEMOD==1){
             newdiam=sqrt((secRa.x(sec)*maxL*4)/(mrgri.x(sec)*PI*100))
             eqdiam.set(index, eqdiam.x(index)+newdiam^2)
        }
    } 
    if(LENGTHMERGINGparallelMETHOD!=3) {eqL.set(index, eqL.x(index)/normalFact)}



    if (PRESERVINGMETHOD==0){             //PRESERVINGMETHOD = 0 DIAMMERGINGparallelMETHOD=0 by Destexhe
       if (DIAMMERGINGparallelMETHOD==-1){
          eqdiam.set(index, sqrt(eqdiamPar.x(index)))
       }else if (DIAMMERGINGparallelMETHOD==0){
          eqdiam.set(index, sqrt(eqdiam.x(index)))
       }else if (DIAMMERGINGparallelMETHOD==1){
          eqdiam.set(index, sqrt(sqrt(eqdiamPar.x(index))*eqdiamSer.x(index)/nsectOrigin))
       }       
    }else if (PRESERVINGMETHOD==1){
       eqdiam.set(index, totsurf/(eqL.x(index)*PI))
    } 

    eqri2.set(index, eqri2Sum/(nsect+1))
    eqSurf.set(index,eqL.x(index)*PI*eqdiam.x(index))
    orSurf1.set(index,surfSum)
    //print "mergingAllMethod: ",totsurf,eqSurf.x(index)

}