      Program doublecounting pulse
      IMPLICIT NONE
      INTEGER N, i, t, Per, gam_t, n_cyc,j,on_part, mem, k, c, z, l  
      INTEGER dist, fc, tc,mc, cont, reacc, noisy,on_noise, cont_n
      INTEGER cont_off, nc, q, reacc_n,on_noise_old,cont_n_old
      INTEGER tot_nbli, loop
      INTEGER state(200000),hist(200000),pho(200000)
      INTEGER dc(200000), n_bli(100),accexp(100000),blink(1000)
      INTEGER bleach(200000),obs(200000),espo(1000,200000), esp(200000)
      INTEGER on_t(200000), tonn(1000),off_t(200000), toffn(1000) 
      real*8   x_tot, y_tot, x_a, y_a, gam_i, lam_i, p_att, ran
      real*8   p_reatt, res, er,check, distx, disty
      real*8   x(200000),y(200000), ranx, rany, checkb, checkn,poff_n 
      real*8   dirt, p_r_n, p_reatt_d,p_off, ransignx, ransigny 

CCCCCCCCCCCCCCC Input parameters  CCCCCCCCCCCCCCC
CCCCCCCCCCCCCCC                   CCCCCCCCCCCCCCC

       
C    Number of particles
       N = 8000
C    volume (x and y) in nm
       x_tot = 3500.0
       y_tot = 3500.0

C    Read the 2D spatial point pattern from an external file

       do i=1,N
        open(unit=10,file='pattern',status='old')
        read(10,'(2F24.14)') x(i), y(i)
       enddo

C    Start cycle over t_d(q)

       do q=1, 100, 1

C    Initialization

       do i=1,200000
        state(i) = 0
        hist(i) = 0
        on_t(i) = 0
        off_t(i) = 0
        pho(i) = 0
        bleach(i) = 0
        dc(i)=0
        obs(i)=0
        esp(i)=0
        do j=1,1000
         espo(j,i)=0
        enddo
       enddo
       do i=1,1000
        tonn(i)=0
        toffn(i)=0
       enddo  
       on_part=0
       tc=0
       mc=0
       cont=0
       reacc=0

CCCCCCCCCCCCCCCCCCC   Experimental settings   CCCCCCCCCCCCCCCCCCC

C    experimental pixel size
       res=130.0
C    data collection memory
       mem=q

C    time (ms)
       t = 0
C    acquistion frame rate (ms)
       Per = 50
c    405 nm illumination cycle 
       loop= 1
c    total number of acquistion frames
       n_cyc= 30000

C    Localization error
      er=10
      

 
CCCCCCCCCCCCCCCCCC Fluorophore photophysics   CCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCC   All these parameters must be tuned on the photophysical behaviour    CCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCC   of the fluorophore in a given environment. The values used in the    CCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCC   following do not represent a particular set-up                       CCCCCCCCCCCCCCCCCCCCC

C    activation probability (related to 405 nm power intensity)

         p_att=0.00015 

C    reactivation probability (tunes t_off)

c       p_reatt= Per/toff 
        p_reatt= 0.11

C    exponential decay Eos (tunes t_on)

c       p_off =  Per/ton 
        p_off = 0.19

C    bleaching probability (depends on t_off, tunes number of blinks) 
        
        checkb = 0.38



CCCCCCCCCCCCCCCCCCC    Noise        CCCCCCCCCCCCCCCCC

C  max number of noise particles that can turn on in a new frame
  
        noisy= 1


C   noise activation probability

        dirt= 0.00 

C   initialization noise state
       
       on_noise=0
       cont_n=0

CCCCCCCCCCCCCCCCCCCCCCCCCCC Start time   CCCCCCCCCCCCCCCCCCCCCCCc

      
      do t=1,n_cyc
       on_part=0

CCCCCCCCCCCCCCCCC Particle dynamics  CCCCCCCCCCCCCCCCCC


C  Switching on  (according to p_att)   
      
      if(mod(t-1,loop)==0) then
       do i=1,N
       if ((state(i).eq.0).and.(hist(i).eq.0).and.(bleach(i).eq.0)) then
       call random_number (ran)
        if (ran.le.p_att) then 
         state(i)=1 
         hist(i)=hist(i)+1
        endif
       endif
       enddo
      endif
 
C  Photoblink         
 
      if(mod(t-1,loop)==0) then 
       do i=1,N    
       if ((state(i).eq.0).and.(hist(i).ne.0).and.(bleach(i).eq.0)) then
        call random_number (ran)
         if (ran.lt.p_reatt) then
          state(i)=1
          hist(i)=hist(i)+1
          toffn(off_t(i))=toffn(off_t(i))+1
          off_t(i)=0
         endif
        endif
        enddo
      else
       do i=1,N
       if ((state(i).eq.0).and.(hist(i).ne.0).and.(bleach(i).eq.0)) then
        call random_number (ran)
         if (ran.lt.p_reatt) then
          state(i)=1
          hist(i)=hist(i)+1
          toffn(off_t(i))=toffn(off_t(i))+1
          off_t(i)=0
         endif
        endif
       enddo

      endif

C  Reversible switching off 

         do i=1,N
           if ((state(i).eq.1)) then
            call random_number (ran)
            if (ran.lt.p_off) then
            state(i) = 0
            tonn(on_t(i))=tonn(on_t(i))+1
            on_t(i) = 0
           endif
          endif
         enddo



C  Particle update

         do i=1,N
          if (state(i).gt.0) then
           on_t(i)=on_t(i)+1
           pho(i)=pho(i)+1
           on_part=on_part+1
          endif
      if ((state(i).eq.0).and.(hist(i).ge.1).and.(bleach(i).eq.0)) then
           off_t(i)=off_t(i)+1
          endif
         enddo


C   Photobleaching 
 
         do i=1,N
          if ((on_t(i).eq.1)) then
          call random_number (ran)
           if (ran.gt.checkb) then
            bleach(i) = 1
           endif
          endif
         enddo

CCCCCCCCCCCCCCCCC End of particle dynamics CCCCCCCCCCCC


CCCCCCCCCCCCCCCCC Start noise dynamics  CCCCCCCCCCCCCCCCCCC


C  Initialization


          on_noise_old=on_noise
          cont_n_old=cont_n  



C  Turn off noise
          cont_off=0
          do i=1,on_noise
            call random_number (ran)
            if (ran.lt.poff_n) then
            cont_off=cont_off+1
            endif
          enddo
           on_noise=on_noise-cont_off



C Turn on noise
 
         do i=1,noisy 
          call random_number (ran)
          if (ran.lt.dirt) then 
           cont_n=cont_n+1
           on_noise=on_noise+1
           call random_number (x_a)
           x(N+cont_n) = x_a*x_tot
           call random_number (y_a)
           y(N+cont_n) = y_a*y_tot
          endif
         enddo


CCCCCCCCCCCCCCCCC End noise dynamics   CCCCCCCCCCCCCCCCCCCCCCCC

CCCCCCCCCCCCCCCCC Experimental (PALM) observation CCCCCCCCCCCCC


C  Initialization

       
         obs(1)=t
         obs(2)=on_part+on_noise
         esp(1)=t
         esp(2)=on_part+on_noise
         c = 1
          do k=1,N
           if (state(k).eq.1) then
            obs(c+2)=k
            c=c+1
           endif
          enddo
          do k=1,on_noise-tot_nbli
           obs(on_part+2+k)=N+cont_n-on_noise+k
          enddo
          do k=1,reacc_n
           if(n_bli(k).gt.0.5) then
           obs(on_part+on_noise-tot_nbli+2+k)=
     &    N+cont_n_old-on_noise_old-k
           endif
          enddo
          open(unit=12,file='test',status='unknown')
        write(12,'(100000I7)') (obs(i), i=1,obs(2)+2)

C  Check, for a given t_d(mem), if the on particles are new localization or belong to old localization       

           do z=1,min(t-1,mem)
           backspace(13)
           enddo
           do z=1,min(t-1,mem)

           read(13,*) (espo(t-(min(t-1,mem))+z-1,i), i=1,
     &                accexp(t-(min(t-1,mem))+z-1)+2) 
           enddo

C  Check whether some of the 'new' localization are due to particles that have already been counted


            do k=3, obs(2)+2
             do z=1,mem
              do l=3, espo(t-z,2)+2
               dist=N
               distx=2*res
               disty=2*res
               distx=abs(int(x(obs(k))/res)-int(x(espo(t-z,l))/res)) 
               disty=abs(int(y(obs(k))/res)-int(y(espo(t-z,l))/res))
               dist=obs(k)-(espo(t-z,l))
               if ((distx.le.1).and.(disty.le.1)) then
                esp(k)=espo(t-z,l)
               go to 300
              endif 
              enddo
             enddo
            esp(k)=obs(k)  
            dc(esp(k))=dc(esp(k))+1

c            write(*,*) t
 
C  Localization output  

           open(unit=11,file='PALM',status='old')
           call random_number (ransignx)
           call random_number (ransigny)
           call random_number (ranx)
           call random_number (rany)
          if(ransignx.gt.0.5) then
           if(ransigny.gt.0.5) then
           write(11,'(2F24.14)') x(esp(k))+ranx*er,y(esp(k))+rany*er
           else
           write(11,'(2F24.14)') x(esp(k))+ranx*er,y(esp(k))-rany*er
           endif
          else
           if(ransigny.gt.0.5) then
           write(11,'(2F24.14)') x(esp(k))-ranx*er,y(esp(k))+rany*er
           else
           write(11,'(2F24.14)') x(esp(k))-ranx*er,y(esp(k))-rany*er
           endif
          endif
 
300         continue
            enddo
        accexp(t)=esp(2)
        open(unit=13,file='esp',status='old')
        write(13,'(100000I7)') (esp(i) ,i=1,esp(2)+2)
       
       enddo
       
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC OUTPUTS  CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


        close(13) 
       nc=0
       do i=1,N+cont_n
        reacc=reacc+hist(i)
        cont=cont+dc(i)
        if(dc(i).gt.1) tc=tc+dc(i)-1
        if(dc(i).lt.0.5) mc=mc+1
        if((i.gt.N).and.(dc(i).gt.0.5)) nc=nc+1
       enddo
       
C      blinks histogram

       do i=1,N
        blink(hist(i))=blink(hist(i))+1
       enddo
       do i=1,1000
        open(unit=25,file='blinks',status='unknown')
        write(25,'(2I9)') i, blink(i)
       enddo 

C      t_on histogram

      do i=1,1000
        open(unit=26,file='t_on',status='unknown')
        write(26,'(2I9)') i*Per, tonn(i)
       enddo

c      t_off histogram

       do i=1,1000
        open(unit=27,file='t_off',status='unknown')
        write(27,'(2I9)') i*Per, toffn(i)
       enddo

C      Summary

       open(unit=14,file='output',status='old')
        write(14,'(6I7)') mem, N,cont, nc, tc, mc 
       enddo

100    end
 
      
