anchor.diff:  
       omegaL(ii)=sqrt(omegae**2.0_wp+3*k(ii)*ve**2.0_wp)
       gamas= (-1)*sqrt(pi/8)*(sqrt(me/mi)+ Te/Ti**2.0_wp/sqrt(Te/Ti)*exp((-1)*(Te/2.0/Ti)-1.5))*abs(k(ii))*Cs
  -		!gamas= (-1)*sqrt(pi/2)*(sqrt(me/mi)+4*pow(Te/2/Ti,2)/sqrt(Te/2/Ti)*exp((-1)*(Te*4/Ti)))*abs(k(ii))*Cs*10;   //based on Robinson 2002
  -		!gamas= (-1)*sqrt(pi/8)*pow(1/(1+k(ii)*k(ii)*lambdaD*lambdaD)+3*Ti/Te,2)/sqrt(1/(1+k(ii)*k(ii)*lambdaD*lambdaD)+3*Ti/Te)*(sqrt(me/mi)+pow(Te/Ti,2)/sqrt(Te/Ti)*exp((-1)*(Te/2.0/Ti)/(1+k(ii)*k(ii)*lambdaD*lambdaD)-1.5))*abs(k(ii))*Cs;   //Based on some Chinese paper!!
  +		!gamas= (-1)*sqrt(pi/2)*(sqrt(me/mi)+4*pow(Te/2/Ti,2)/sqrt(Te/2/Ti)*exp((-1)*(Te*4/Ti)))*abs(k(ii))*Cs*10   //based on Robinson 2002
  +		!gamas= (-1)*sqrt(pi/8)*pow(1/(1+k(ii)*k(ii)*lambdaD*lambdaD)+3*Ti/Te,2)/sqrt(1/(1+k(ii)*k(ii)*lambdaD*lambdaD)+3*Ti/Te)*(sqrt(me/mi)+pow(Te/Ti,2)/sqrt(Te/Ti)*exp((-1)*(Te/2.0/Ti)/(1+k(ii)*k(ii)*lambdaD*lambdaD)-1.5))*abs(k(ii))*Cs   //Based on some Chinese paper!!
       nui(ii)=(nuic/2-gamas)
   
       if (ii==N/2) then
  
positive.code: 
          do q=LL,UU
            CC(1)=CC(1)+(EE(c1,q+N/2,1)+k2(q+N/2,1)/2)*(EE(c1,q-p(pp)+N/2,1)+k2(q-p(pp)+N/2,1)/2) + &
                  (EE(c1,q+N/2,2)+k2(q+N/2,2)/2)*(EE(c1,q-p(pp)+N/2,2)+k2(q-p(pp)+N/2,2)/2);
            CC(2)=CC(2)+(EE(c1,q+N/2,2)+k2(q+N/2,2)/2)*(EE(c1,q-p(pp)+N/2,1)+k2(q-p(pp)+N/2,1)/2) - &
                  (EE(c1,q+N/2,1)+k2(q+N/2,1)/2)*(EE(c1,q-p(pp)+N/2,2)+k2(q-p(pp)+N/2,2)/2);
          end do
  
          kn3(:)=Tstep*(vv(c1,pp,:)+kv2(:)/2+SSn(:)/2*Tstep)
  
negative.code: 
        do pp=1,N/2
          call random_number(rdist(:2))
          SSn(:)= rdist(:2)*Source_factor_n(pp)/sqrt(Tstep)
  
          LL= max(p(pp)-N/3,-N/3)
          UU= min(N/3,p(pp)+N/3)
          CC(:)=0.0_wp
  
          do q=LL,UU
            ! no vect
            CC(1)=CC(1)+EE(c1,q+N/2,1)*EE(c1,q-p(pp)+N/2,1)+EE(c1,q+N/2,2)*EE(c1,q-p(pp)+N/2,2)
            CC(2)=CC(2)+EE(c1,q+N/2,2)*EE(c1,q-p(pp)+N/2,1)-EE(c1,q+N/2,1)*EE(c1,q-p(pp)+N/2,2)
          end do
  
          kn1(:)=Tstep*(vv(c1,pp,:))
  