diff --git a/src/starformation.c b/src/starformation.c
index d1ae4c2..779b24d 100644
--- a/src/starformation.c
+++ b/src/starformation.c
@@ -1,1236 +1,1238 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
 #include <mpi.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_spline.h>
 
 #include "allvars.h"
 #include "proto.h"
 
 #ifdef PY_INTERFACE
 #include <Python.h>
 #include <numpy/arrayobject.h>
 #endif
 
 #ifdef SFR
 #ifdef FEEDBACK_WIND 
 
 void feedbackwind_compute_energy_kin(int mode)
 {
   int i;
   double DeltaEgyKin;
   double Tot_DeltaEgyKin;
 
   DeltaEgyKin = 0;
   Tot_DeltaEgyKin = 0;
   
   if (mode==1)
     {
       LocalSysState.EnergyKin1 = 0;
       LocalSysState.EnergyKin2 = 0;
    }
     
     
   for(i = 0; i < N_gas; i++)
     {
       if (P[i].Type==0)
         {
     
       	  if (mode==1)
             LocalSysState.EnergyKin1 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]);
       	  else
             LocalSysState.EnergyKin2 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]);
         }
     }	
     
   if (mode==2)   
     {
       DeltaEgyKin = LocalSysState.EnergyKin2 - LocalSysState.EnergyKin1;
       MPI_Reduce(&DeltaEgyKin, &Tot_DeltaEgyKin, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); 
       LocalSysState.EnergyFeedbackWind -= DeltaEgyKin;
     }
  
       
 }
 
 #endif
 
 
 
 
 
 
 /*! \file starformation.c
  *  \brief Compute jet
  *
 */
 
 static double hubble_a,a3inv=1.;
 
 /*! Compute starformation
  *
  */
 void star_formation(void)
 {
   int i,j;
   int NumNewStars=0,firststar;
   int *AllNumNewStars;
   double T;
   double p;
   double dt=0;
   double tstar=1.;
   double MeanBaryonicDensity;
   int flagST=1,flagToStar=0;
   double mstar;
   double mstars=0,Tot_mstars=0;
   double mnewstars,Tot_mnewstars=0;
   int Nnewstars,Tot_Nnewstars=0;
   double star_formation_rate;
   int *numlist;
   double t0,t1;
   double DivVel;
   
   struct particle_data Psave;
   
   int k;
 
   t0 = second();
 
   if(ThisTask == 0)
     {
       printf("start star formation... \n");
       fflush(stdout);
     }  
     
   AllNumNewStars = malloc(sizeof(int) * NTask);
   
 
 #ifdef FEEDBACK_WIND  
   double phi,costh,sinth;
 #endif  
 
   if(All.ComovingIntegrationOn)
     {
       hubble_a = All.Omega0 / (All.Time * All.Time * All.Time)
 	+ (1 - All.Omega0 - All.OmegaLambda) / (All.Time * All.Time) + All.OmegaLambda;
 
       hubble_a = All.Hubble * sqrt(hubble_a);
       a3inv = 1 / (All.Time * All.Time * All.Time);
     }
   else
     {
       hubble_a = 1.;
       a3inv = 1.;
     }
         
   switch (All.StarFormationType)  
     {
       case 0:
         All.ThresholdDensity = All.StarFormationDensity;
         break;
 	
       case 1:
         MeanBaryonicDensity = All.OmegaBaryon * (3 * All.Hubble * All.Hubble / (8 * M_PI * All.G));
         All.ThresholdDensity = dmax(All.StarFormationDensity/a3inv,200*MeanBaryonicDensity);  
         break;	
 	
       case 2:
         All.ThresholdDensity = All.StarFormationDensity;
         All.StarFormationTime = 1./All.StarFormationCstar / pow(4*PI*All.G*All.StarFormationDensity,0.5);
         break;	
         
       case 3:
         All.StarFormationTime = 1./All.StarFormationCstar / pow(4*PI*All.G*All.StarFormationDensity,0.5);
         /* All.ThresholdDensity is fixed further */
         break;          
     }
   
   flagST=1;
   mnewstars=0;
   Nnewstars=0;
 #ifdef FEEDBACK_WIND  
   feedbackwind_compute_energy_kin(1);
 #endif
     
   for(i = 0; i < N_gas; i++)	       
     {     	
       /* only active particles */  
       if(P[i].Ti_endstep == All.Ti_Current)
         {
           
 #ifdef CHIMIE_KINETIC_FEEDBACK	
         if (SphP[i].WindTime < (All.Time-All.ChimieWindTime))	/* not a wind particle */
 	  {
 #endif
 
 	  /* if the particle is a false gas (new star) */
 	  if(P[i].Type != 0)
 	    continue;
 
 #ifdef FEEDBACK		
 	  /* if the particle is a SN */	  
           if(SphP[i].EnergySN)
 	    { 
 	      
 	      /* the particle is not allowed to becomes a stars, as it is already a SN */
 	      flagST = 0;	
 	      
 	      SphP[i].EnergySNrem = SphP[i].EnergySNrem - SphP[i].EnergySN;
 	     
 	      /* dt in time unit here, even for ComovingIntegrationOn */
               dt = (All.Ti_Current - P[i].Ti_begstep) * All.Timebase_interval / hubble_a;
 
 	      if (All.SupernovaTime==0)
                 SphP[i].EnergySN = SphP[i].EnergySNrem;
 	      else 
                 SphP[i].EnergySN = (P[i].Mass*All.SupernovaEgySpecPerMassUnit)/All.SupernovaTime*dt;	
 	      
 	      
 	      /* limit the feedback energy to the remaining energy */
               if (SphP[i].EnergySN>SphP[i].EnergySNrem)
                 SphP[i].EnergySN = SphP[i].EnergySNrem; 	    
 	      
 	      
 	      if (SphP[i].EnergySNrem == 0)		/* not enough energy, make a star */
 	        {
 		  //printf("Particle %d becomes a star\n",P[i].ID);
 		  P[i].Type = ST;
 		  RearrangeParticlesFlag=1;
                 }
 
 		
 	    }
 #endif	    
 		
 	  /* check if the particle may become a star */ 					        
 #ifdef SFR_NEG_DIV_ONLY
           if(All.ComovingIntegrationOn)
 	    DivVel = SphP[i].DivVel/(All.Time*All.Time) + 3*hubble_a;
 	  else
 	    DivVel = SphP[i].DivVel;
 	  
 	  if ((  DivVel   < 0) && flagST) 	
 #endif	  							        
 	    {		
 	    	    									        
 
          
 	 
 	 
 #ifdef MULTIPHASE	 
 	      if (SphP[i].Phase == GAS_STICKY)
 	        {
 #endif
                   if (All.StarFormationType==3)
                     All.ThresholdDensity = 0.5* ( pow(SphP[i].DivVel,2) + pow(SphP[i].CurlVel,2)  )/(All.G);
 
 	      	  if (SphP[i].Density*a3inv > All.ThresholdDensity)
 	      	    {	
 		    	
                                             
 	      	      /* check temperature */								      								      
 #ifdef ISOTHERM_EQS  
               	      T = All.InitGasTemp;
 #else         	  		
 #ifdef MULTIPHASE
                       T = GAMMA_MINUS1*All.mumh/All.Boltzmann * SphP[i].EntropyPred;
 #else
                       T = All.mumh/All.Boltzmann * SphP[i].EntropyPred * pow(SphP[i].Density*a3inv,GAMMA_MINUS1); 
 #endif
 #endif	      	  											      
 
 	      	      if (T < All.StarFormationTemperature)								      
 	      	  	{	    
 		  	  /* dt in time unit here, even for ComovingIntegrationOn			 */
               	  	  /* dt has units of Gyr/h also tstar						 */
 	      	  	  dt = (All.Ti_Current - P[i].Ti_begstep) * All.Timebase_interval / hubble_a;	    
 		  	  tstar = All.StarFormationTime / pow(SphP[i].Density*a3inv/All.StarFormationDensity,0.5);
 		  	  
 			  
 			  /* check the mass of the particle */
 			  
 			  flagToStar = 0;
 			  
 			  if (All.StarFormationNStarsFromGas==1)
 			    {
 			      /* transform the gas particle into a star particle */ 
 			      mstar = P[i].Mass;
 			      flagToStar=1;
 			    }
 			  else
 			    {  
 			      if (P[i].Mass<=(1.0+All.StarFormationMgMsFraction)*All.StarFormationStarMass)
 			        /* transform the gas particle into a star particle */  
 			        {
 			           mstar = P[i].Mass;
 				   flagToStar=1;
 			        }
                               else
 			        /* extract a star particle from the gas particle */
 			        {
 			           mstar = All.StarFormationStarMass;
 				   flagToStar=0;
 			        }
 			    }
 			    
 			    
 			    
 			  
 			  /* compute probability of forming star */
 	      	  	  p = (P[i].Mass/mstar)*(1. - exp(-dt/tstar) ); 
                           		  														       
       	      	  	  if (get_StarFormation_random_number(P[i].ID) < p)
 		  	    {
                                                                                                                                               
                               			      			    	    
 #ifdef FEEDBACK
                               SphP[i].EnergySNrem = P[i].Mass*All.SupernovaEgySpecPerMassUnit;
 			      
 			      if (All.SupernovaTime==0)
                                 SphP[i].EnergySN    = SphP[i].EnergySNrem;
 		              else
                                 SphP[i].EnergySN    = (P[i].Mass*All.SupernovaEgySpecPerMassUnit)/All.SupernovaTime*dt;
 			      
 			      			      
 			      /* limit the feedback energy to the remaining energy */
 			      if (SphP[i].EnergySN>SphP[i].EnergySNrem)
 			        SphP[i].EnergySN = SphP[i].EnergySNrem; 
 						              			      		  	      
 		  	      if (SphP[i].EnergySNrem==0)
 		  		{
 		  		  P[i].Type = ST;
 		  		  RearrangeParticlesFlag=1;
 		  		}
 #else	
 
                               if (flagToStar)	/* turns the stellar particule into star particle */
                                 {
 #ifdef STELLAR_PROP
 				  if (N_stars+1 > All.MaxPartStars)
    				    {
    				      printf("No memory left for new star particle : N_stars+1=%d MaxPartStars=%d\n\n",N_stars+1,All.MaxPartStars);
    				      endrun(999008);
    				    } 	
 #else
 				  if (N_stars+1 > All.MaxPart-N_gas)
    				    {
    				      printf("No memory left for new star particle : N_stars+1=%d All.MaxPart-N_gas=%d\n\n",N_stars+1,All.MaxPart-N_gas);
    				      endrun(999009);
    				    } 	
 
 #endif				    
                                   	           						                  
 				  P[i].Type = ST;	
 			          RearrangeParticlesFlag=1;
 				  
 				  //printf("(%d) move gas particle to star (i=%d,id=%d)\n",ThisTask,i,P[i].ID);
 				  
 				  /* count mass */
 				  mnewstars+=P[i].Mass;
 				  Nnewstars+=1;
 				  
 				  				  				  				  
 				  /* gather energy */
 #ifdef FEEDBACK
 	      			  LocalSysState.StarEnergyFeedback += SphP[i].EgySpecFeedback * P[i].Mass;
 #endif 				  
 				  
 				}	  
 			      else	/* create a new star particle */
 			        { 
 				  			  
 #ifdef STELLAR_PROP				  
 				  if (N_stars+1 > All.MaxPartStars)
    				    {
    				      printf("No memory left for new star particle : N_stars+1=%d MaxPartStars=%d\n\n",N_stars+1,All.MaxPartStars);
    				      endrun(999000);
    				    } 	
 #else
 				  if (N_stars+1 > All.MaxPart-N_gas)
    				    {
    				      printf("No memory left for new star particle : N_stars+1=%d All.MaxPart-N_gas=%d\n\n",N_stars+1,All.MaxPart-N_gas);
    				      endrun(999001);
    				    } 	
 #endif				    			  
 				  				  
 				  if (NumPart+1 > All.MaxPart)
    				    {
    				      printf("No memory left for new particle : NumPart=%d MaxPart=%d",NumPart,All.MaxPart);
    				      endrun(999002);
    				    } 
 				    	   				  
 				  if (P[i].Mass <= mstar)
 				    {
    				      printf("Mass too small %g %g !",P[i].Mass,mstar);
    				      endrun(999003);				    
 				    }
 				 
    				 
    				  /* the new particle get all properties of the SPH particle, except the mass */
    				  /* the type is defined further */
    				 
    				  j = NumPart+NumNewStars;			  
    				 
    				  P[j]         = P[i];
    				  P[j].Mass    = mstar;
    				  P[j].Type    = ST;
 
 				  /* count mass */
 				  mnewstars+=mstar;
 				  Nnewstars+=1;
 
 				     				 
    				  /* this is needed to kick parent nodes dynamically in timestep.c */
    				  /* could be avoided if All.NumForcesSinceLastDomainDecomp > All.TotNumPart * All.TreeDomainUpdateFrequency*/
    				  Father[j] = Father[i];
 				     				 
 #ifdef STELLAR_PROP				
 				  /* index to Stp */
                                   P[j].StPIdx = N_stars+NumNewStars;
 
 				  /* record info at the end of StP */
   
 				  /* record time */
 				  StP[P[j].StPIdx].FormationTime = All.Time;
 				  /* record time */
 				  StP[P[j].StPIdx].IDProj        = P[i].ID;   
 				  /* record hsml */    
 				  StP[P[j].StPIdx].Hsml          = SphP[i].Hsml; 
 
 #ifdef CHIMIE
 				  /* record initial mass */    
 				  StP[P[j].StPIdx].InitialMass   = P[j].Mass; 
 				  /* record density */    
 				  StP[P[j].StPIdx].Density       = SphP[i].Density; 	
 				  /* record metalicity */    
 	        		  for(k = 0; k < NELEMENTS; k++)
                         	    StP[P[j].StPIdx].Metal[k]    = SphP[i].Metal[k];
 #endif 				  			  
 				  /* index to P */
 				  StP[P[j].StPIdx].PIdx		 = j;
                                    
 				  
 #endif	
 
 
 				  /* change proj properties */
    				  P[i].Mass = P[i].Mass-mstar;
 
 
 			  
 				  
 				  /* gather energy */
 #ifdef FEEDBACK
 	      			  LocalSysState.StarEnergyFeedback += SphP[i].EgySpecFeedback * P[j].Mass;
 #endif 
 
 				  /* here, we should add the self force, but it is negligible */
 
 
    				  NumNewStars++;
 
 				 				 
 			       }	 
 			      			      			       
 #endif				
 		  	    } 
 #ifdef FEEDBACK_WIND			    
 			  else
 			    {
 			      tstar = tstar/All.SupernovaWindParameter;
 			      /* here, we devide by p, because there is (1-p) prob. to be here */
 			      p = (P[i].Mass/mstar)*(1. - exp(-dt/tstar))/(1-p);	
 			      
 			      if (get_FeedbackWind_random_number(P[i].ID) < p)
 			        {
                                   phi	= PI*2.*get_FeedbackWind_random_number(P[i].ID+1);	
 		  	   	  costh = 1.-2.*get_FeedbackWind_random_number(P[i].ID+2);	
 		  		  sinth = sqrt(1.-costh*costh);    
 		  
 				  SphP[i].FeedbackWindVel[0] = All.SupernovaWindSpeed  *sinth*cos(phi);	
 		   		  SphP[i].FeedbackWindVel[1] = All.SupernovaWindSpeed  *sinth*sin(phi);	
 		   		  SphP[i].FeedbackWindVel[2] = All.SupernovaWindSpeed  *costh; 
 				  
 				  
 				  //printf("(%d) %d vx vy vz %g %g %g   %g\n",ThisTask,P[i].ID,P[i].Vel[0],P[i].Vel[1],P[i].Vel[2],sqrt( P[i].Vel[0]*P[i].Vel[0] + P[i].Vel[1]*P[i].Vel[1] + P[i].Vel[2]*P[i].Vel[2] ));
 				  //printf("(%d) SupernovaWindSpeed   %g   dv=%g\n",ThisTask,All.SupernovaWindSpeed,sqrt( SphP[i].FeedbackWindVel[0]*SphP[i].FeedbackWindVel[0] + SphP[i].FeedbackWindVel[1]*SphP[i].FeedbackWindVel[1] + SphP[i].FeedbackWindVel[2]*SphP[i].FeedbackWindVel[2] ));
 				  
 			 	  /* new velocities */
 	  	      	 	  for(k = 0; k < 3; k++)
       	  	      	 	    {				   
                                       SphP[i].VelPred[k] += SphP[i].FeedbackWindVel[k];
                                       P[i].Vel[k]        += SphP[i].FeedbackWindVel[k];
 			 	    }	
 				    
 				  //printf("(%d) %d vx vy vz %g %g %g   %g\n",ThisTask,P[i].ID,P[i].Vel[0],P[i].Vel[1],P[i].Vel[2],sqrt( P[i].Vel[0]*P[i].Vel[0] + P[i].Vel[1]*P[i].Vel[1] + P[i].Vel[2]*P[i].Vel[2] ));
 				  
 				  
 		                }  
 			 	
 				
 				
 			    }
 #endif			    
 			      
 			     
 
 	      	  	}									      
               	    }
 
 #ifdef MULTIPHASE		
 		}
 #endif			
 		
 		
 		
 													  
 	    }											        
 
 
 #ifdef CHIMIE_KINETIC_FEEDBACK	
 	  }
 #endif
 
 
 	}
 
     }
     
 
 
 
   if (All.StarFormationNStarsFromGas != 1)
     {
      
           
       /* get NumNewStars from all other procs */
       MPI_Allgather(&NumNewStars, 1, MPI_INT, AllNumNewStars, 1, MPI_INT, MPI_COMM_WORLD);
           
       /* find the id of the fist local star */
       firststar =  All.MaxID; 
   
       for (i=0;i<ThisTask;i++)
         firststar += AllNumNewStars[i];
   
       /* now, set ids */  
       for (i=0;i< NumNewStars;i++)
         P[NumPart+i].ID = firststar + i; 
 	              
       /* move the new stars to the right place
         in both cases, the relati ve position of the new stars is preseved, 
         so, there is no need to rearrange StP.
       */
       if (NumNewStars < NumPart-N_gas-N_stars)			/* move stars inwards */
         {
 	  for(i=0;i<NumNewStars;i++)
 	    {	    
               Psave = P[N_gas+N_stars+i];     
 	      P[N_gas+N_stars+i] = P[NumPart+i];
               P[NumPart+i]=Psave;	
 #ifdef STELLAR_PROP   	      
 	      /* set index for StP */
 	      StP[P[N_gas+N_stars+i].StPIdx].PIdx = N_gas+N_stars+i; 
 #ifdef CHECK_ID_CORRESPONDENCE	
               StP[P[N_gas+N_stars+i].StPIdx].ID = P[N_gas+N_stars+i].ID; 
 #endif   
               //printf("(%d) new star (a) i=%d id=%d (StPi=%d PIdx=%d)\n",ThisTask,N_gas+N_stars+i,P[N_gas+N_stars+i].ID,P[N_gas+N_stars+i].StPIdx,StP[P[N_gas+N_stars+i].StPIdx].PIdx);
 #endif 	      
 
 	    }
 	} 
       else							/* move other particles outwards */
         { 
 	
 	  if (NumPart-N_gas-N_stars>0)				/* there is other particles */
 	    {
 	      for(i=N_gas+N_stars;i<NumPart;i++)
 	        {		
                   Psave = P[i]; 	  
 	          P[i] = P[i+NumNewStars];
 	          P[i+NumNewStars] = Psave;
 #ifdef STELLAR_PROP 	  
 	          /* set new index for StP */
 	          StP[P[i].StPIdx].PIdx = i; 
 #ifdef CHECK_ID_CORRESPONDENCE  
                   StP[P[i].StPIdx].ID = P[i].ID; 
 #endif          
 	          
 #endif 	 
                   //printf("(%d) new star (b) i=%d id=%d\n",ThisTask,N_gas+N_stars+i,P[N_gas+N_stars+i].ID);
 	        }
 
 
 
 	      /* now, set new index for the StP that where untouched */	
 	      for(i=N_gas+N_stars+1;i<N_gas+N_stars+NumNewStars;i++)		/* here, we should move up to <N_gas+N_stars+NumNewStars */
 	        {	                                                        /* however, if there is only one particle, it will not be counted */
 #ifdef STELLAR_PROP 	  
 	          /* set new index for StP */
 	          StP[P[i].StPIdx].PIdx = i; 
 #ifdef CHECK_ID_CORRESPONDENCE  
                   StP[P[i].StPIdx].ID = P[i].ID; 
 #endif          
 #endif 	
                 }
 
 
 
 
             }
 	  else							/* there is no other particles outside, i.e, new particles are already the last ones */
 	    {
 	    
 	      for(i=N_gas+N_stars;i<NumPart+NumNewStars;i++)
 	        {	    
 #ifdef STELLAR_PROP 	  
 	          /* set new index for StP */
 	          StP[P[i].StPIdx].PIdx = i; 
 #ifdef CHECK_ID_CORRESPONDENCE  
                   StP[P[i].StPIdx].ID = P[i].ID; 
 #endif          
 	          
 #endif 	 
 	        }
 	    
 	      
 	    }
 	    
 	    
 
 
 
+
         }
 
 
 
 
       /* set new NumPart */
       NumPart = NumPart+NumNewStars;
       N_stars = N_stars+NumNewStars;
 
 
 
       /*MPI_Allreduce(&NumPart,  &All.TotNumPart,  1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);*/
       numlist = malloc(NTask * sizeof(int) * NTask);
 
       MPI_Allgather(&NumPart, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD); 
       for(i = 0, All.TotNumPart = 0; i < NTask; i++)
         All.TotNumPart += numlist[i];
 	
       MPI_Allgather(&N_stars, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD); 
       for(i = 0, All.TotN_stars = 0; i < NTask; i++)
         All.TotN_stars += numlist[i];
 
       free(numlist);  
       
+      
+      All.MaxID += AllNumNewStars;      
 
     }
 
 
   /* compute star mass */
   mstars=0;
   for (i=0;i< NumPart;i++)
     {
       if (P[i].Type==ST)
         mstars+= P[i].Mass;
     }
     
   
 #ifdef FEEDBACK_WIND  
   feedbackwind_compute_energy_kin(2);
 #endif
 
 
 
   /* share results */
   MPI_Allreduce(&mstars,  &Tot_mstars,  1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
   MPI_Allreduce(&mnewstars,  &Tot_mnewstars,  1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
   MPI_Allreduce(&Nnewstars,  &Tot_Nnewstars,  1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
   
-  All.MaxID += Tot_Nnewstars; 
   
   if (All.TimeStep>0)
     star_formation_rate = Tot_mnewstars/All.TimeStep;
   else
     star_formation_rate = 0.;
 
 #ifndef PY_INTERFACE
   if (ThisTask==0)
     {  
       //fprintf(FdSfr, "Step %d, Time: %g, NewStarsPart: %g \n", All.NumCurrentTiStep, All.Time, Tot_mstars);
       fprintf(FdSfr, "%15g %15g %15g %15g %15g\n", All.Time,All.TimeStep,Tot_mstars,Tot_mnewstars,star_formation_rate);
       fflush(FdSfr);
     }  
 #endif
 
  //printf("(%d) N_gas = %d  N_stars = %d NumPart = %d\n",ThisTask,N_gas,N_stars,NumPart);
  //fflush(stdout);
 
  if(ThisTask == 0)
    { 
      printf("Total number of new star particles :  %d    \n",(int) Tot_Nnewstars);
      printf("Total number of star particles     :  %d%09d\n",(int) (All.TotN_stars / 1000000000), (int) (All.TotN_stars % 1000000000));
      fflush(stdout);
    }
 
 
 
 #ifdef CHECK_BLOCK_ORDER
   int typenow;
 
   rearrange_particle_sequence();
 
   typenow = 0;
   for(i = 0; i < NumPart; i++)    
     {
     	       
       if ((P[i].Type<typenow)&&(typenow<=ST))
     	{
     	  printf("\nBlock order error\n");
     	  printf("(%d) i=%d id=%d Type=%d typenow=%d\n\n",ThisTask,i,P[i].ID,P[i].Type,typenow);
     	  endrun(999004);
     	}
     		 
      typenow = P[i].Type;     
     }	    
 #endif
 
 
 
 #ifdef CHECK_ID_CORRESPONDENCE
 
   rearrange_particle_sequence();
 
   for(i = N_gas; i < N_gas+N_stars; i++)    
     {
 
       if( StP[P[i].StPIdx].PIdx != i )
        	{
     	  printf("\nP/StP correspondance error\n");
     	  printf("(%d) (in starformation) N_stars=%d N_gas=%d i=%d id=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d\n\n",ThisTask,N_stars,N_gas,i,P[i].ID,P[i].StPIdx,StP[P[i].StPIdx].PIdx);
     	  endrun(999005);
     	}	        
       
 
       if(StP[P[i].StPIdx].ID     != P[i].ID)
        	{
     	  printf("\nP/StP correspondance error\n");
     	  printf("(%d) (in starformation) N_gas=%d N_stars=%d i=%d Type=%d P.Id=%d P[i].StPIdx=%d StP[P[i].StPIdx].ID=%d \n\n",ThisTask,N_gas,N_stars,i,P[i].Type,P[i].ID, P[i].StPIdx, StP[P[i].StPIdx].ID);
     	  endrun(999006);
     	}		
 	       
     }	    
 
 #endif
 
   free(AllNumNewStars);
 
 
   if(ThisTask == 0)
     {
       printf("star formation done. \n");
       fflush(stdout);
     }  
     
   
   t1 = second();  
   All.CPU_StarFormation += timediff(t0, t1);
   
 }
 
 
 /*
  * This routine rearrange particle sequence in order to group particles of same type
  */
 void rearrange_particle_sequence(void)
 {
   int i,j;
   struct particle_data Psave;
   struct sph_particle_data Psave_sph;  
   int nstars=0;
   int *numlist;
   int sumflag;
 #ifdef CHIMIE
   int k;
 #endif
   
 
   MPI_Allreduce(&RearrangeParticlesFlag, &sumflag, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); 
 
  
 
   if(RearrangeParticlesFlag)
     {
        
     
       if(ThisTask == 0)
         {
           printf("Start rearrange particle sequence...\n");
           fflush(stdout);
         }    
     
     
       /* loop over all gasous particles */
       for(i = 0; i < N_gas; i++)
         {	
 	
           if (P[i].Type == ST)
             {      
 
 	      
 	      /* find the first non star (flag) particle */
 	      for(j=N_gas-1;j>0;j--)
 	        {
                   if (P[j].Type==0)	/* we have a gasous particle */
 		    break;
 
 		  if (j==i)		/* the particule is the last gaseous particle */
 		    break;		/* the particle will be inverted with itself */
 
 				  
 		  if (P[j].Type==ST) 	/* the particule is selected to becomes a star, turn it to star */
 		    {
 #ifdef STELLAR_PROP
 		      P[j].StPIdx = N_stars;
 		    
               	      StP[P[j].StPIdx].FormationTime = All.Time;  		
               	      StP[P[j].StPIdx].IDProj 	     = P[j].ID;
 	      	      StP[P[j].StPIdx].Hsml          = SphP[j].Hsml;
 #ifdef CHIMIE
                       /* record initial mass */    
                       StP[P[j].StPIdx].InitialMass   = P[j].Mass; 
 		      /* record density */    
 		      StP[P[j].StPIdx].Density       = SphP[j].Density;     
 		      /* record metalicity */	 
 	              for(k = 0; k < NELEMENTS; k++)
                         StP[P[j].StPIdx].Metal[k]    = SphP[j].Metal[k];		     
 #endif 				  			  
 
 #ifdef CHECK_ID_CORRESPONDENCE	
                       StP[P[j].StPIdx].ID            = P[j].ID; 
 #endif  
 		      	      
 		      StP[P[j].StPIdx].PIdx          = j;
                       //printf("(%d) move to stars (b) i=%d -> j=%d id=%d Stpi=%d (N_gas=%d N_stars=%d)\n",ThisTask,j,j,P[j].ID,P[j].StPIdx,N_gas,N_stars);
 #endif		      
 
               	      nstars++;
               	      N_gas--;     
 	      	      N_stars++;  
 		      continue;   			     
 		    }
 
 
 		}
 	      
 	      //printf("(%d) --> i=%d id=%d Type=%d (N_gas=%d)\n",ThisTask,i,P[i].ID,P[i].Type,N_gas);
 	      //printf("(%d) <-- i=%d id=%d Type=%d\n",ThisTask,j,P[j].ID,P[j].Type);
 
 	      /* move the particle */	    
               Psave = P[i];									      	/* copy the particle */
 	      Psave_sph = SphP[i];								      	/* copy the particle */
 	      
 	      P[i] = P[j];										/* replace by the last gaseous one */
 	      SphP[i] = SphP[j];									/* replace by the last gaseous one */								
 	      P[j] = Psave;									
 	      
 #ifdef STELLAR_PROP
               P[j].StPIdx = N_stars;
               StP[P[j].StPIdx].FormationTime = All.Time;				
               StP[P[j].StPIdx].IDProj        = P[j].ID;
 	      StP[P[j].StPIdx].Hsml          = Psave_sph.Hsml;
 #ifdef CHIMIE
               /* record initial mass */    
               StP[P[j].StPIdx].InitialMass   = P[j].Mass; 
 	      /* record density */    
 	      StP[P[j].StPIdx].Density       = Psave_sph.Density;     
 	      /* record metalicity */	 
 	      for(k = 0; k < NELEMENTS; k++)
                 StP[P[j].StPIdx].Metal[k]    = Psave_sph.Metal[k];		     
 #endif 	
 
 
 #ifdef CHECK_ID_CORRESPONDENCE	
               StP[P[j].StPIdx].ID            = P[j].ID;
 #endif  
 	      StP[P[j].StPIdx].PIdx	     = j;
               //printf("(%d) move to stars (a) i=%d -> j=%d id=%d Stpi=%d (N_gas=%d N_stars=%d)\n",ThisTask,i,j,P[j].ID,P[j].StPIdx,N_gas,N_stars);
 #endif
 
 
 
               nstars++;
               N_gas--;	   
 	      N_stars++;   	  	      
             }
         }
       RearrangeParticlesFlag=0;
     
       //if(ThisTask == 0)
       //  {
       //    printf("%d new star(s)\n",nstars);
       //    printf("Start rearrange particle sequence done.\n");
       //    fflush(stdout);
       //  }    
     
 
     }
     
    
 
     if (sumflag)		/* do this only if at least one Task has rearrange particles */
       {
       
       	
 
       	numlist = malloc(NTask * sizeof(int) * NTask);
 
       	/* update all.TotN_gas */
       	MPI_Allgather(&N_gas, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD); 
       	for(i = 0, All.TotN_gas = 0; i < NTask; i++)
       	  All.TotN_gas += numlist[i];
 
       	/* update all.TotN_stars */
       	MPI_Allgather(&N_stars, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD); 
       	for(i = 0, All.TotN_stars = 0; i < NTask; i++)
       	  All.TotN_stars += numlist[i];
 
       	free(numlist);  
 
 
       }
 
 
 
 
 
 }
 
 
 
 
 
 
 void sfr_compute_energy_int(int mode)
 {
   int i;
   double DeltaEgyInt;
   double Tot_DeltaEgyInt;
   double egyspec;
 
 
   DeltaEgyInt = 0;
   Tot_DeltaEgyInt = 0;
   
   if (mode==1)
     {
       LocalSysState.EnergyInt1 = 0;
       LocalSysState.EnergyInt2 = 0;
    }
         
   for(i = 0; i < N_gas; i++)
     {
       if (P[i].Type==0)
         {
 
 #ifdef MULTIPHASE
           if (SphP[i].Phase== GAS_SPH)
 	    egyspec = SphP[i].Entropy / (GAMMA_MINUS1) * pow(SphP[i].Density * a3inv, GAMMA_MINUS1);
 	  else
 	    egyspec = SphP[i].Entropy;
 #else
 	  egyspec = SphP[i].Entropy / (GAMMA_MINUS1) * pow(SphP[i].Density * a3inv, GAMMA_MINUS1);
 
 #endif
     
       	  if (mode==1)
             LocalSysState.EnergyInt1 += P[i].Mass * egyspec;
       	  else
             LocalSysState.EnergyInt2 += P[i].Mass * egyspec;
         }
     }	
     
   if (mode==2)   
     {
       DeltaEgyInt = LocalSysState.EnergyInt2 - LocalSysState.EnergyInt1;
       MPI_Reduce(&DeltaEgyInt, &Tot_DeltaEgyInt, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD); 
       LocalSysState.StarEnergyInt -=  DeltaEgyInt;   
     }
  
       
 }
 
 
 
 
 
 
 
 
 
 
 
 
 
 void sfr_check_number_of_stars(int mode)
 {
 
   int i;
   int nstars;
   int npotstars;
   
          
   for(i = 0, nstars = 0, npotstars = 0; i < NumPart; i++)
     {
        
       if (P[i].Type==ST)
         if (i<N_gas)
 	  {
 	    npotstars++;
 	  }
 	else
 	  {
 	    nstars++;
 	  } 
 	
     }
 
 
   printf("(%d) N_stars=%d  npotstars+nstars=%d npotstars=%d nstars=%d\n",ThisTask,N_stars,npotstars+nstars,npotstars,nstars);
   if (nstars != N_stars)
     endrun(987987);
 
 }
 
 
 
 /****************************************************************************************/
 /*
 /*
 /*
 /*                               PYTHON INTERFACE                                     
 /*
 /*
 /*
 /****************************************************************************************/
 
 
 
 /*********************************/
 /*                               */
 /*********************************/
       
 #ifdef PY_INTERFACE
 
 PyObject * sfr_SetParameters(PyObject *self, PyObject *args)
 {
         
     PyObject *dict;
     PyObject *key;
     PyObject *value;
     int ivalue;
     float fvalue;
     double dvalue;
             
         
     /* here, we can have either arguments or dict directly */
     if(PyDict_Check(args))
       {
         dict = args;
       }
     else
       {
         if (! PyArg_ParseTuple(args, "O",&dict))
           return NULL;
       } 
         
     /* check that it is a PyDictObject */
     if(!PyDict_Check(dict))
       {
         PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary.");   
         return NULL;
       }
     
      Py_ssize_t pos=0;
      while(PyDict_Next(dict,&pos,&key,&value))
        {
          
          if(PyString_Check(key))
            {
                    
                                  
              if(strcmp(PyString_AsString(key), "SfrFile")==0)
                {
                  if(PyString_Check(value))
                    strcpy(  All.SfrFile , PyString_AsString(value));                                   
                }   
                
                
              if(strcmp(PyString_AsString(key), "StarsAllocFactor")==0)
                {
                  if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
                      All.StarsAllocFactor = PyFloat_AsDouble(value);                                   
                }                 
              
              if(strcmp(PyString_AsString(key), "StarFormationNStarsFromGas")==0)
                {
                  if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
                      All.StarFormationNStarsFromGas = PyInt_AsLong(value);                                   
                }        
                
              if(strcmp(PyString_AsString(key), "StarFormationMgMsFraction")==0)
                {
                  if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
                      All.StarFormationMgMsFraction = PyFloat_AsDouble(value);                                   
                }                   
 
                     
              if(strcmp(PyString_AsString(key), "StarFormationStarMass")==0)
                {
                  if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
                      All.StarFormationStarMass = PyFloat_AsDouble(value);                                   
                } 
  
              if(strcmp(PyString_AsString(key), "StarFormationType")==0)
                {
                  if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
                      All.StarFormationType = PyInt_AsLong(value);                                   
                } 
 
               if(strcmp(PyString_AsString(key), "StarFormationCstar")==0)
                {
                  if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
                      All.StarFormationCstar = PyFloat_AsDouble(value);        
                }    
 
               if(strcmp(PyString_AsString(key), "StarFormationTime")==0)
                {
                  if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
                      All.StarFormationTime = PyFloat_AsDouble(value);                                   
                }       
     
               if(strcmp(PyString_AsString(key), "StarFormationDensity")==0)
                {
                  if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
                      All.StarFormationDensity = PyFloat_AsDouble(value);                                   
                }   
 
               if(strcmp(PyString_AsString(key), "StarFormationTemperature")==0)
                {
                  if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
                      All.StarFormationTemperature = PyFloat_AsDouble(value);                                   
                }   
 
                                                
            }
        }
 
 
     return Py_BuildValue("i",1);
 }
 
 
 
 
 PyObject * sfr_GetParameters()
 {
   
     PyObject *dict;
     PyObject *key;
     PyObject *value; 
     
     dict = PyDict_New();   
     
     key   = PyString_FromString("SfrFile");
     value = PyString_FromString(All.SfrFile);
     PyDict_SetItem(dict,key,value);
 
     key   = PyString_FromString("StarsAllocFactor");
     value = PyFloat_FromDouble(All.StarsAllocFactor);
     PyDict_SetItem(dict,key,value);   
             
     key   = PyString_FromString("StarFormationNStarsFromGas");
     value = PyInt_FromLong(All.StarFormationNStarsFromGas);
     PyDict_SetItem(dict,key,value);  
 
     key   = PyString_FromString("StarFormationMgMsFraction");
     value = PyFloat_FromDouble(All.StarFormationMgMsFraction);
     PyDict_SetItem(dict,key,value);  
     
     key   = PyString_FromString("StarFormationStarMass");
     value = PyFloat_FromDouble(All.StarFormationStarMass);
     PyDict_SetItem(dict,key,value);     
     
     key   = PyString_FromString("StarFormationType");
     value = PyInt_FromLong(All.StarFormationType);
     PyDict_SetItem(dict,key,value);     
     
     key   = PyString_FromString("StarFormationCstar");
     value = PyFloat_FromDouble(All.StarFormationCstar);
     PyDict_SetItem(dict,key,value);      
     
     key   = PyString_FromString("StarFormationTime");
     value = PyFloat_FromDouble(All.StarFormationTime);
     PyDict_SetItem(dict,key,value);      
     
     key   = PyString_FromString("StarFormationDensity");
     value = PyFloat_FromDouble(All.StarFormationDensity);
     PyDict_SetItem(dict,key,value);    
   
     key   = PyString_FromString("StarFormationTemperature");
     value = PyFloat_FromDouble(All.StarFormationTemperature);
     PyDict_SetItem(dict,key,value);  
     
     return Py_BuildValue("O",dict);
     
 }
 
 
 
 
 
 
       
 
 
 
 
 /* definition of the method table */      
       
 static PyMethodDef sfrMethods[] = {
   
   
           {"SetParameters",  sfr_SetParameters, METH_VARARGS,
            "set parameters"},    
 
           {"GetParameters",  sfr_GetParameters, METH_VARARGS,
            "get parameters"},  
   
 
           {NULL, NULL, 0, NULL}        /* Sentinel */
       };      
       
 
 
 
 void initsfr(void)
       {    
           (void) Py_InitModule("sfr", sfrMethods);  
           
           import_array();
       }      
 
 
 
 
 
 #endif /* PY_INTERFACE */
 #endif /* SFR */