polsim1.txt

This file contains the documentation and source code for the 
Fortran program referred to in Montgomery & Swank, ApJ 801,21
and included in the Technical Report available from
http://asd.gsfc.nasa.gov/Jean.Swank/

The program generates and analyses data sets like those in X-ray
polarimetry measurements.  It is written in Standard Fortran77, 
and so is also Standard F90 and F95 and F2003 and F2008, and 
should work with no problems on any contemporary fortran system.  
It is short and simple, and should be easy to convert, if desired, 
to another language -- IDL, Matlab, Python, Haskell, ...

The random number generator used is the RAND function, which is 
not standard but almost universally provided by Fortran 
compilers. It could, and probably should, be replaced by the 
standard F95 subroutines RANDOM_SEED and RANDOM_NUMBER, which are 
also available on any fortran system less than 20 years old. 
Replacement with any other random number generator would be 
simple. This application does not need a particularly good 
generator; any default is probably adequate.

The function IPOISSON is a minimal one, adequate for this use. 
There are many more sophisticated ones available. If the number 
of counts per bin gets up into the hundreds an improvement might 
be desirable (or the number of bins increased). The limitation 
comes from underflow in the exponential of a large negative 
number, and is therefore quite system-dependent.

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c  polsim1.f  
c  simulation of counts and analysis for x-ray polarization 
c  cgm  2013 

c  Standard Fortran77 except for "implicit none" and "rand" 
c
c  Variable names intended to resemble those in text 
c  Nbar is intended value for number of counts 
c  Ntot is actual value in a particular sample 
c  Phiest is angle estimate in degrees 
c  Uz, Qz are displacements from point Z 
c  R2 is square of distance from Z 
c  D2 is square of ellipse semi-major axis 
c 
c  Usage: choose values for Nsamp,Nbar,Azero 
c         adjust write statement if desired 
c         compile, link, execute.
c

      implicit none
      integer Mbins, Nsamp, Ntot, ipoisson
      parameter (Mbins = 499)
      integer i, k, ks 
      integer nj(Mbins), seed(3)
      real Nbar, Azero, OneoverPi, p, f, g, R
      real Vcos(Mbins), Vsin(Mbins)
      real Ic, Uc, Qc, Aest, Phiest
      real Uz, Qz, R2, D2, lam2

      Nsamp = 49000
      Nbar = 4000.
      Azero = 0.50

      OneoverPi = 0.25/atan(1.0)
      ks = (Mbins+1)/2
      f = 16*atan(1.0)/Mbins
      g = Nbar/Mbins
      lam2 = 1.0/(1.0-0.5*Azero**2)

      DO 10 k=1,Mbins
        p = f*(k-ks)
        Vcos(k) = cos(p)
        Vsin(k) = sin(p)
  10  CONTINUE

      CALL itime(seed)
      R=rand(seed(1)*seed(2)+seed(3))
      OPEN (13, file='psimout')

      DO 50 i=1,Nsamp
        ntot = 0
        DO 20 k=1,Mbins
          R = g*( 1.0+Azero*Vcos(k) )
          nj(k) = ipoisson(R)
          Ntot = Ntot + nj(k)
  20    CONTINUE

        Uc = 0.0
        Qc = 0.0
        DO 30 k=1,Mbins
          Uc=Uc + nj(k)*Vcos(k)
          Qc=Qc + nj(k)*Vsin(k)
  30    CONTINUE

        Uc=OneoverPi * Uc
        Qc=OneoverPi * Qc
        Ic = 0.5*OneoverPi*Ntot
        Uz=Uc/Ic-Azero
        Qz=Qc/Ic
        R2=Uz**2 + Qz**2
        D2=Qz**2 + lam2*Uz**2
        Aest = sqrt(Uc**2 + Qc**2)/Ic
        Phiest = 90*OneoverPi*atan2(Qc,Uc)

        WRITE(13,*) Uz, Qz, Aest, Phiest, Ic, Uc, Qc

  50  CONTINUE

      CLOSE(13)
      END

C********************************

       function ipoisson(a) 
C----------returns integer from Poisson distribution with mean a 
C------------uses intrinsic rand(flag), which is probably C srand 
C------------can be slow for arguments much larger than 10
C------------not reliable for arguments over 400 or so

      real a
      double precision b, c
      integer i
      b=exp(-1.d0*a)
      i=-1
      c=1
   10 i=i+1
      c=c*rand(0)
      if (c .gt. b) goto 10
      ipoisson=i
      return
      end 

cccccccc end of file ccccccccccccccccc