text: import numpy
  
  def ZernikeAnnularJacobian(Z, x, y , e, atype):
  
  
      m1, n1 = x.shape
      m2, n2 = y.shape
  
      if( m1 != m2 or n1 != n2 ):
          print( 'x & y are not the same size' )
          exit()
  
      if( len(Z) > 22 ):
          print('ZernikeAnnularEval() is not implemented with >22 terms')
          return
      elif len(Z)<22:
          Z[21]=0
  
      x2 = x* x
      y2 = y* y
      xy = x* y
      r2 = x2 + y2
      x4=x2*x2
      x6=x4*x2
      y4=y2*y2
      y6=y4*y2
      e2=e*e
      e4=e2*e2
      e6=e4*e2
      e8=e6*e2
      e10=e8*e2
      e12=e10*e2
      e14=e12*e2
      e16=e14*e2
      
      if (atype == '1st'):
          j =        Z[0]  * 0 * x # to make d an array with the same size as x
          j =  j + Z[1] * 0
          j =  j + Z[2]  * 0
          den=1-e**2
          j =  j + Z[3]  * numpy.sqrt( 3 )  * 8/den
          j =  j + Z[4] * numpy.sqrt( 6 )  * 0
          j =  j + Z[5] * numpy.sqrt( 6 )  * 0
          den=numpy.sqrt( (1-e2)**2*(1+e2)*(1+4*e2+e4) )
          j =  j + Z[6] * numpy.sqrt( 8 )  * 24 * y*(1+e2)/den 
          j =  j + Z[7] * numpy.sqrt( 8 )  * 24 * x*(1+e2)/den 
          j =  j + Z[8] * numpy.sqrt( 8 )  * 0
          j =  j + Z[9] * numpy.sqrt( 8 )  * 0
          den=(1-e2)**2
          j =  j + Z[10] * numpy.sqrt( 5 )  * (96 * r2 - 24*(1+e2))/den
          den=(1-e2)**3*(1+e2+e4)
          num=numpy.sqrt((1-e2)**4*(1+e2+e4)/(1+4*e2+10*e4+4*e6+e8))
          j =  j + Z[11] * numpy.sqrt( 10 ) * 48 * (x2 - y2)*(1+e2+e4)*num/den 
          j =  j + Z[12] * numpy.sqrt( 10 ) * 96 * xy*(1+e2+e4)*num/den
          j =  j + Z[13] * numpy.sqrt( 10 ) * 0
          j =  j + Z[14] * numpy.sqrt( 10 ) * 0
          den=(1-e2)**3*(1+4*e2+e4)
          num=numpy.sqrt((1-e2)**2*(1+4*e2+e4)/(1+9*e2+9*e4+e6) )
          j =  j + Z[15] * numpy.sqrt( 12 ) * 48*x*(5*r2*(1+4*e2+e4)-2*(1+4*e2+4*e4+e6))*num/den
          j =  j + Z[16] * numpy.sqrt( 12 ) * 48*y*(5*r2*(1+4*e2+e4)-2*(1+4*e2+4*e4+e6))*num/den
          den=(1-e2)**4*(1+e2)*(1+e4)
          num=numpy.sqrt((1-e2)**6*(1+e2)*(1+e4)/(1+4*e2+10*e4+20*e6+10*e8+4*e10+e12))    
          j =  j + Z[17] * numpy.sqrt( 12 ) * 80.0*x*(x2-3.0*y2)*(1+e2)*(1+e4)*num/den
          j =  j + Z[18] * numpy.sqrt( 12 ) * 80.0*y*(3*x2-y2)*(1+e2)*(1+e4)*num/den
          j =  j + Z[19] * numpy.sqrt( 12 ) * 0
          j =  j + Z[20] * numpy.sqrt( 12 ) * 0
          den=(1-e2)**3
          j =  j + Z[21] * numpy.sqrt(  7 )  * 48*(e4 -10*e2*x2 - 10*e2*y2 + 3*e2 + 15*x4 + 30*x2*y2- 10*x2 +15*y4-10*y2+1)/den
      elif (atype =='2nd'):
      
          j =        Z[0]**2  * 0 * x # to make d an array with the same size as x
          j =  j + Z[1]**2  * 0
          j =  j + Z[2]**2  * 0
          den=1-e**2
          j =  j + Z[3]**2 * ( 3 )  * 16/den/den
          den=(1+e2+e4)    
          j =  j + Z[4]**2 * ( 6 )  * (-4)/den
          j =  j + Z[5]**2 * ( 6 )  * (-4)/den
          den=(1-e2)**2*(1+e2)*(1+4*e2+e4)
          j =  j + Z[6]**2 * ( 8 )  * (108 * y2 - 36 * x2)*(1+e2)/den 
          j =  j + Z[7]**2 * ( 8 )  * (108 * x2 - 36 * y2)*(1+e2)/den 
          den=(1+e2+e4+e6)
          j =  j + Z[8]**2 * ( 8 )  * (-36 * r2)/den
          j =  j + Z[9]**2 * ( 8 )  * (-36 * r2)/den
          den=(1-e2)**4
          j =  j + Z[10]**2 * ( 5 )  * 144 * (1+e2-2*r2)*(1+e2-6*r2)/den
          den=(1-e2)**6*(1+e2+e4)**2
          num=((1-e2)**4*(1+e2+e4)/(1+4*e2+10*e4+4*e6+e8))
          j =  j + Z[11]**2 * ( 10 ) * 36 * (8 *(1+e2+e4)* x2 - 1-e2-e4-e6)* (1+e2+e4+e6 - 8 *(1+e2+e4)* y2)*num/den
          j =  j + Z[12]**2 * ( 10 ) * 36 * (-4*(x-y)**2*(e4+e2+1)+1+e2+e4+e6)*(4*(x+y)**2*(e4+e2+1)-1-e2-e4-e6)*num/den
          den=(1+e2+e4+e6+e8)
          j =  j + Z[13]**2 * ( 10 ) * (-144) * r2**2/den
          j =  j + Z[14]**2 * ( 10 ) * (-144) * r2**2/den
          den=(1-e2)**6*(1+4*e2+e4)**2
          num=(1-e2)**2*(1+4*e2+e4)/(1+9*e2+9*e4+e6)
          j =  j + Z[15]**2 * ( 12 ) * 64*((3*e6 - 5*e4*r2 + 12*e4 - 20*e2*r2 + 12*e2 - 5*r2  + 3)*\
                                           (9*e6*x2 - 3*e6*y2 - 25*e4*x4 - 20*e4*x2*y2 + 36*e4*x2 + 5*e4*y4 - 12*e4*y2 - 100*e2*x4 - 80*e2*x2*y2\
                                            + 36*e2*x2 + 20*e2*y4 - 12*e2*y2 - 25*x4 - 20*x2*y2 + 9*x2 + 5*y4 - 3*y2)  )*num/den
          j =  j + Z[16]**2 * ( 12 ) * 64*(-(3*e6 - 5*e4*r2 + 12*e4 - 20*e2*r2 + 12*e2 - 5*r2 + 3)* \
                                            (3*e6*x2 - 9*e6*y2 - 5*e4*x4 + 20*e4*x2*y2 + 12*e4*x2 + 25*e4*y4 - 36*e4*y2 - 20*e2*x4 + 80*e2*x2*y2\
                                             + 12*e2*x2 + 100*e2*y4 - 36*e2*y2 - 5*x4 + 20*x2*y2 + 3*x2 + 25*y4 - 9*y2)  )*num/den
          den=(1-e2)**8*(1+e2)**2*(1+e4)**2
          num=(1-e2)**6*(1+e2)*(1+e4)/(1+4*e2+10*e4+20*e6+10*e8+4*e10+e12)
          j =  j + Z[17]**2 * ( 12 ) * 16.0*(- 36*e16*x2 - 36*e16*y2 + 180*e14*x4 + 360*e14*x2*y2 - 72*e14*x2\
                                               + 180*e14*y4 - 72*e14*y2 - 125*e12*x6 - 1275*e12*x4*y2 + 360*e12*x4 + 225*e12*x2*y4\
                                               + 720*e12*x2*y2 - 108*e12*x2 - 225*e12*y6 + 360*e12*y4 - 108*e12*y2 - 250*e10*x6 \
                                               - 2550*e10*x4*y2 + 540*e10*x4 + 450*e10*x2*y4 + 1080*e10*x2*y2 - 144*e10*x2 - 450*e10*y6\
                                               + 540*e10*y4 - 144*e10*y2 - 375*e8*x6 - 3825*e8*x4*y2 + 720*e8*x4 + 675*e8*x2*y4\
                                               + 1440*e8*x2*y2 - 180*e8*x2 - 675*e8*y6 + 720*e8*y4 - 180*e8*y2 - 500*e6*x6 - 5100*e6*x4*y2\
                                               + 720*e6*x4 + 900*e6*x2*y4 + 1440*e6*x2*y2 - 144*e6*x2 - 900*e6*y6 + 720*e6*y4 - 144*e6*y2\
                                               - 375*e4*x6 - 3825*e4*x4*y2 + 540*e4*x4 + 675*e4*x2*y4 + 1080*e4*x2*y2 - 108*e4*x2 - 675*e4*y6\
                                               + 540*e4*y4 - 108*e4*y2 - 250*e2*x6 - 2550*e2*x4*y2 + 360*e2*x4 + 450*e2*x2*y4 + 720*e2*x2*y2\
                                               - 72*e2*x2 - 450*e2*y6 + 360*e2*y4 - 72*e2*y2 - 125*x6 - 1275*x4*y2 + 180*x4 + 225*x2*y4\
                                               + 360*x2*y2 - 36*x2 - 225*y6 + 180*y4 - 36*y2 )*num/den
          j =  j + Z[18]**2 * ( 12 ) * 16.0*((- 225*e12 - 450*e10 - 675*e8 - 900*e6 - 675*e4 - 450*e2 - 225)*x6\
                                             + (180*e14 + 225*e12*y2 + 360*e12 + 450*e10*y2 + 540*e10 + 675*e8*y2 + 720*e8 + 900*e6*y2\
                                                + 720*e6 + 675*e4*y2 + 540*e4 + 450*e2*y2 + 360*e2 + 225*y2 + 180)*x4\
                                             + (- 36*e16 + 360*e14*y2 - 72*e14 - 1275*e12*y4 + 720*e12*y2 - 108*e12 - 2550*e10*y4\
                                                  + 1080*e10*y2 - 144*e10 - 3825*e8*y4 + 1440*e8*y2 - 180*e8 - 5100*e6*y4 + 1440*e6*y2\
                                                  - 144*e6 - 3825*e4*y4 + 1080*e4*y2 - 108*e4 - 2550*e2*y4 + 720*e2*y2 - 72*e2 - 1275*y4\
                                                  + 360*y2 - 36)*x2 - 36*e16*y2 + 180*e14*y4 - 72*e14*y2 - 125*e12*y6 + 360*e12*y4 - 108*e12*y2\
                                             - 250*e10*y6 + 540*e10*y4 - 144*e10*y2 - 375*e8*y6 + 720*e8*y4 - 180*e8*y2 - 500*e6*y6\
                                             + 720*e6*y4 - 144*e6*y2 - 375*e4*y6 + 540*e4*y4 - 108*e4*y2 - 250*e2*y6 + 360*e2*y4 - 72*e2*y2\
                                             - 125*y6 + 180*y4 - 36*y2 )*num/den
          den=(1+e2+e4+e6+e8+e10)
          j =  j + Z[19]**2 * ( 12 ) * (-400) * r2**3/den
          j =  j + Z[20]**2 * ( 12 ) * (-400) * r2**3/den
          den=(1-e2)**6
          j =  j + Z[21]**2 * (  7 )  * 576*((e4 - 5*e2*x2 - 5*e2*y2 + 3*e2 + 5*x4 + 10*x2*y2 - 5*x2 + 5*y4 - 5*y2 + 1)* \
                                             (e4 - 15*e2*x2 - 15*e2*y2 + 3*e2 + 25*x4 + 50*x2*y2 - 15*x2 + 25*y4 - 15*y2 + 1) )/den
  
  
      return j
  