In this part we defined the relevant constants.  alpha is where the pulse begins and bet where it ends in the 
theta coordinates.  Period is the period of the cycle in the new coordinates, i.e. the periodicity of phi.

alph  =  0.18*Pi;  bet  =  0.28*Pi;  l  =  1.736; 
    DD  = 1.05;

Next we define the period of the cycle and the value of the phi coordinate when it ends.  Note that in phi coordinates
the perturbation lasts from alph to bdry which makes phi T periodic rather than 2*Pi periodic.

period = 2*Pi  +  l/(1 + l) * (alph  - bet)

bdry  =  (bet  +  l*alph)/(1 + l);

The following definitions then give phi as a function of theta.  The first is valid between alph and bet, while the second is
valid after bet.  This just follows from the condition phi' = te

def1[s_]  :=  (1 + l)*s  -  alph*l;  
    def2[s_]  :=  s  +  (bet - alph) * l/(1 + l);

Next we define the function corresponding to sin 2*theta in the paper. This is a function of phi, so that it is different in the interval
[alph,bdry] and outside of it.  We also take into account the different definitions in each of the separate regions.

p[t_ ]  :=  (s  =  Mod[t,  period]; 
        Which[s > alph  &&  
            s <  bdry,  -Sin[2*def1[s]]/(l + 
                1), 
                                                                             s  >  
            bdry,  -Sin[
              2*def2[s]],  
                                                                             True,   -Sin[
              2*s]])

Plot[p[t],  {t,  0,  period}]

Next we define the places where the singularities in the integration occur.  This is necessary because otherwise the integrator
will complain.

n = 0;  sing =  {t, 0};
alphmod  =  (alph  -  Psi)/DD;
 bdrymod  =  (bdry  - Psi)/DD;
         While[(n*period/DD + alphmod) < 
        period/DD,  
                  If[(n*period/DD + alphmod) >  0 ,   sing = 
          Append[sing, 
            n*period/DD + 
              alphmod]]; 
                  If [(n*period/DD + bdrymod) >  
            0  &&  (n*period/DD + bdrymod)  <  
            period/DD,                sing = 
          Append[sing, 
            n*period/DD + 
              bdrymod]
    ];  
  n  =  n + 1
]
sing = Append[sing,  period/DD];

Lastly we define the integral.  This procedure just computes the points of discontinuity first, puts them in a list, and then integrates over them.
The list sing contains the values where the singularities appear.

V[Psi_]  :=  (n = 0;  
        sing =  {t, 
            0}; 
alphmod  =  (alph  -  Psi)/
            DD; 
 bdrymod  =  (bdry  - Psi)/
                DD;
         While[(n*
                    period/DD + alphmod) < 
              period/DD,  
                  If[(n*period/DD + alphmod) >  0 ,   sing = 
                Append[sing, 
                  n*period/DD + 
                    alphmod]]; 
                  If [(n*period/DD + bdrymod) >  
                  0  &&  (n*period/DD + bdrymod)  <  
                  period/DD,                sing = 
                Append[sing, 
                  n*period/DD + 
                    bdrymod]
    ];  
  n  =  
              n + 1
]
sing = 
            Append[sing,  period/DD];  
 
        1/period*  NIntegrate[p[DD*t + Psi], 
            Evaluate[sing]] )

Now we can compute the contribution from the term sin 2*theta on the dynamics of the phase.   The fluctuations are due to numerical integration.
This could have been done in a much more straigthforward fashion, but it allowed us to set up the functions
necessary in the next section.

Plot[0.075*V[Psi]/(2*Pi),  {Psi,  0,  2*Pi}]

Since this term is small and constant for all values of Psi, it will only shift the Arnold tongue slightly.

Calculating the Gamma function

Next we define the omega function and integrate as above.  R denotes the radius, while o is the driving frequency.

R =  11;  o  =  (2*Pi) / (DD*period);

Here we compute the approximation to the function Gamma, by ignoring the jumps.

od/(2*Pi)*Integrate[-Sin[od t + Psi] * 
          Sin[od*t]/R,  {t,  0,  2*Pi/od}]

s2  =  Plot[- Cos[t]/(11*2),  {t,  0,  2*Pi}]

Next we compute the actual value of the function by taking into account the jumps.  Again, we compute the set of 
singularities.   The function Omega is defined differently on three intervals.

Omega[phi_]  :=  (s  =  Mod[phi,  period]; 
        Which[s > alph  &&  s <  bdry,  -Sin[def1[s]]*
            Sin[o*t] /(R*(l + 
                      1)), 
                                                                             s  >  
            bdry,  -Sin[def2[s]] *
            Sin[o*t]/
              R,  
                                                                             True,   -Sin[s]*
            Sin[o*t]/R])

gamma[Psi_,  DD_]  :=  (n = 0;  
        sing =  {t, 0}; 
period = 
          2*Pi  +  
            l/(1 + l) * (alph  - 
                  bet); 
bdry  =  (bet  +  
                l*alph)/(1 + 
                l); 
alphmod  =  (alph  -  Psi)/
            DD; 
 bdrymod  =  (bdry  - Psi)/
                DD;
         While[(n*
                    period/DD + alphmod) < 
              period/DD,  
                  If[(n*period/DD + alphmod) >  0 ,   sing = 
                Append[sing, 
                  n*period/DD + 
                    alphmod]]; 
                  If [(n*period/DD + bdrymod) >  
                  0  &&  (n*period/DD + bdrymod)  <  
                  period/DD,                sing = 
                Append[sing, 
                  n*period/DD + 
                    bdrymod]
    ];  
  n  =  
              n + 1
]
sing = 
            Append[sing,  period/DD];  
   (1/
              period)*
          NIntegrate[Omega[DD*t + Psi], Evaluate[sing]] )

Now we plot the exact solution in red, and the approximate solution from above in black and compare the two.

s1 =  
      Plot[gamma[t,  1],  {t,  0,  period}, 
        PlotStyle -> {RGBColor[1, 0, 0]}]

Show[s1, s2]

Next we compute the Arnold tongue.  First we find the points at which the Gamma reaches its minimum or maximum depending on the sign of epsilon Delta.

r =.

min[DD_] :=  
      Mod[If [DD < 1,  
          r  /. FindMinimum[
                gamma[r,  DD],  {r,  1.5,  
                  5}][[2]],  
r  /. FindMinimum[-gamma[r,  DD],  {r,  1.5,  5}][[2]]],  
        period]

The function epsilon now gives the boundary of the Arnold tongue.

epsilon[DD_]  :=  (DD - 1)/gamma[min[DD],  DD]

Finally we compare the results with the data.

p2  = 
      Plot[epsilon[DD + 1],  {DD ,  -0.004,  0.005}]

Next we define the data that was calculated from numerical experiments:

data  =  ReadList["nb12.dat",  {Number,  Number}]

data1 = 
        Table[{data[[n, 1]] /1.0329  -  1, data[[n, 2]]}, {n, 1, 
            Length[data]}];

p1  =  
      ListPlot[data1,  Prolog -> AbsolutePointSize[6]]

fig = Show[p1,  p2,  
        Ticks  ->  {Automatic,  {0.025,  0.05,  0.075,  0.1}}, 
        AxesLabel  -> "coupling strength",  
        Epilog  ->   Text["frequency ratio",  {1.003,  
              0.01}]]

Export["arnoldtongue.eps",  
      fig,  "EPS"]