use GD;
package scale;
# This package is used for storing scaling information when plotting using GD. It is not used when generating gcode. 
# Its simple and fairly versatile, but still results in inverted y axis which are not taken account 
# of, so produces upside down plots.  This doesn matter. Much. :-(

$VERSION=0.01; 

sub new
{
    my $s={};

    $s->{ox}=$s->{oy}=0; 
    $s->{s}=1.0; 

    return bless $s;
}

sub setorigin
{
    my ($s,$x,$y)=@_; 

    $s->{ox}=$x;
    $s->{oy}=$y;

}
sub setpixelorigin
{
    my ($s,$x,$y)=@_; 

    $s->{pox}=$x;
    $s->{poy}=$y;

}


sub setscale
{
    my ($s,$scale,$xsize,$ysize)=@_; 

    

    $s->{s}=$xsize/$scale; 
    $s->{xsize}=$xsize; 
    $s->{ysize}=$ysize; 
    $s->{pox}=$xsize/2; 
    $s->{poy}=$ysize/2; 
}

sub scalexy
{
    my $s=shift(@_); 
    my (@xyus)=@_;  # unscaled xy pairs
    my @xys;        # scaled xy pairs
#   print "scalexy b4 is @xyus\n"; 
    while (@xyus)
    {
        my $x=$xyus[0]; 
        my $y=$xyus[1];
        push(@xys, (int(0.5+($x-$s->{ox})*$s->{s})+$s->{pox}),int(0.5+(($y-$s->{oy})*$s->{s})+$s->{poy})); 
        shift(@xyus); shift(@xyus); 
    }
#    print "scalexy after is @xys\n"; 
    return (@xys); 
}
sub scaled # scale a distance only, no origin offset, use for things like diameters
{
    my $s=shift(@_); 
    return map { int(0.5+$_*$s->{s}) } @_;  
}

# This holds the information for a a single wheel: a pinion or spur gear.  
package wheel;
$VERSION=0.01; 
my $pi=3.1415926; 
sub new
{
    my ($s,$m,$n)=@_; 
    $s={};
    $s->{n}=$n; 
    $s->{m}=$m; 

    return bless $s;
}
# Trepan: - to cut a hole in. If you want holes in your wheels, call this functions. ( Ie you want spokes!) 
sub trepan
{
  my ($w,        # pointer to self, a wheel. 
        $spoken, # number of spokes
        $wos,    # total width of spokes in inches
        $bsf,    # boss radius as a factor of pitch radius (dimentionless) 
        $rsf,    # rim size factor as proportion of pitch radius. (dimentionless) 
        $roe,     # radius of window edge   (in inches, the curved radius of the join between a spoke and the rim or boss
                 # not that this must be more than the cutter size or you cant cut it! This is radius, not diameter. 
        $wobf,    # width at base factor, > 1 for wider spoke base. Tapered spokes anyone ? 
        $srf     # spoke rotation factor rotates spoke position this proportion of a full rotation
                  # on outer rim 
                 # values 0 to 0.1 give good results. 
        )=@_; 

# ($spoken,$wos,$bsf,$rsf,$roe,$wobf)=(6,0.5,0.25,0.35,0.075,1.0); 

# for now we just store the values in the wheel . 
        $w->{spoken}=$spoken; 
        $w->{wos}=$wos;  
        $w->{bsf}=$bsf;   
        $w->{rsf}=$rsf;   
        $w->{roe}=$roe;   
        $w->{wobf}=$wobf;    
        $w->{srf}=$srf; 
}
# actually cut the requested trepanning schee

# note that the nasty smooth algorithmn only works if the circle is centered on origin 0,0. 
# for the moment, easy way to correct this is to do all calculations assuming origin based wheel 
# then offset just before plotting with xi yi, the initial position, which is the real wheel center. 
sub cuttrepan
{
    my ($cp,
        $gp, $xi,$yi ) =@_;  

    my ($spoken, # number of spokes
        $wos,    # total width of spokes
        $bsf,    # boss radius as a factor of pitch radius
        $rsf,    # rim size facter as proportion of pitch radius. 
        $roe,     # radius of window edge
        $wobf,   # width at base factor, > 1 for wider spoke base
        $srf
        )= map { $cp->{$_} } qw(spoken wos bsf rsf roe wobf srf); 

    map{ eval "\$$_=$cp->{$_}" } qw(spoken wos bsf rsf roe wobf); 

    return if (!$spoken);   # no spokes, no trepanning. 

    my $mm=$cp->{mm}; 
    my $pi=3.1415926; 
    my $bossradius=$mm*$bsf*$cp->{dw}/2;

    my $rr=$mm*(1-$rsf)*$cp->{dw}/2; # rim radius; 
    my $mm=$cp->{mm}; 

    $wos=0.5*$wos/$spoken;
    $wosb=$wobf*$wos; # width at boss
    
    $gp->gmove('z',0.1); 
    $srf=$srf*2*$pi; 
    my @xy,@l2;
    my $tx,$ty; # temory x,y variables; 

    for $w  ( 0..$spoken-1)   # for each window, we calculate all the points, and put them on a stack. 
    {                         # before we plot, we process to radius the sharp edges. 
        my $t1=2*$pi*$w/$spoken; 
        my $t2=$t1+2*$pi/$spoken; # end of this window
        my $d; 
        $first=1;
        my $z=0; 
        $d=$wosb/$bossradius;

        my $passno=0; 
        while ($passno++ < $cp->{passes})
        {

        if ($first)
        {
          $x=($bossradius+$cp->{cuttersize}/2)*cos($t1+$pi/$spoken);                      # positioning for 2nd half of circle segment
          $y=($bossradius+$cp->{cuttersize}/2)*sin($t1+$pi/$spoken);
          $gp->gmove('x',$xi+$x,'y',$yi+$y);  
        }
        $x=$bossradius*cos($t1+$pi/$spoken);                      # positioning for 2nd half of circle segment
        $y=$bossradius*sin($t1+$pi/$spoken); 
        push(@xy,$x,$y); 
        
        $x+=$bossradius*(cos($t2-$d)-cos($t1+$pi/$spoken));       # rotation at boss radius
        $y+=$bossradius*(sin($t2-$d)-sin($t1+$pi/$spoken));
        my $t=$t2-$d; 
        push(@xy,$x,$y); 
         
        $d=$wos/$rr; 
        $x+=$rr*cos($t2-$d+$srf)-$bossradius*cos($t);                  # radial move to rim 
        $y+=$rr*sin($t2-$d+$srf)-$bossradius*sin($t);
        push(@xy,$x,$y); 
 
        $x+=$rr*(cos($t1+$d+$srf)-cos($t2-$d+$srf));                        # rotatin around rim 
        $y+=$rr*(sin($t1+$d+$srf)-sin($t2-$d+$srf));
        $t=$t1+$d+$srf; 
        push(@xy,$x,$y); 
 
        $d=$wosb/$bossradius; 
        $x+=$bossradius*cos($t1+$d)-$rr*cos($t);                   # radial move back towards center 
        $y+=$bossradius*sin($t1+$d)-$rr*sin($t); 
        push(@xy,$x,$y);

        $x=$bossradius*cos($t1+$pi/$spoken);                       # remaining half of rotation around boss. 
        $y=$bossradius*sin($t1+$pi/$spoken); 
        push(@xy,$x,$y); 

        $z+=$cp->{passdepth}; # passdeth is negative
        
        if ($first)
        {
          $gp->gcompr('d',1,$gp->gmove('x',$xi+shift(@xy),'y',$yi+shift(@xy),'z',$z));     # initial position
                                                                                         # ccomp on, move to z=0; 
        }
        else 
        {
           $gp->gmove('x',$xi+shift(@xy),'y',$yi+shift(@xy),'z',$z)
           # shift(@xy); shift(@xy); 
        }
        $z+=$cp->{passdepth} ; # passdepth -ve 

        @l2=rsmooth(shift(@xy),shift(@xy),shift(@xy),shift(@xy),$bossradius,-$roe);
        @xy=(@l2,@xy); 
        $gp->garcccw('x',$tx=$xi+shift(@xy),'y',$ty=$yi+shift(@xy),'r',$bossradius,'z',$z); # rotation at boss radius, add in z incremen 


        @l2=smooth(shift(@xy),shift(@xy),shift(@xy),shift(@xy),$rr,$roe); 
        @xy=(@l2,@xy);                                                   
        $gp->garccw('x',$xi+shift(@xy),'y',$yi+shift(@xy),'r',$roe);         
       
        $gp->gmove('x',$xi+shift(@xy),'y',$yi+shift(@xy));               # line outwards
        $gp->garccw('x',$xi+shift(@xy),'y',$yi+shift(@xy),'r',$roe);            
       
        
        @l2=rsmooth(shift(@xy),shift(@xy),shift(@xy),shift(@xy),$rr,-$roe);
        @xy=(@l2,@xy); 
        $gp->garccw('x',$xi+shift(@xy),'y',$yi+shift(@xy),'r',$rr);      # outer radius

        @l2=smooth(shift(@xy),shift(@xy),shift(@xy),shift(@xy),$bossradius,$roe); 
        @xy=(@l2,@xy); 

        $gp->garccw('x',$xi+shift(@xy),'y',$yi+shift(@xy),'r',$roe);
        #         
        $gp->gmove('x',$xi+shift(@xy),'y',$yi+shift(@xy));              # line inwards
        $gp->garccw('x',$xi+shift(@xy),'y',$yi+shift(@xy),'r',$roe);
        $gp->garcccw('x',$xi+shift(@xy),'y',$yi+shift(@xy),'r',$bossradius);
#       $gp->gmove('z',0.1); 
        $first=0; 
      } # all passes complete. 
     # repeat 1st move as z was being ramped up during this move. 
     $gp->garcccw('x',$tx,'y',$ty,'r',$bossradius,'z',$z); # rotation at boss radius, add in z incremen  
     $gp->gmove('z',0.1); 
    }
}
sub cuthole
{
  my ($cp,$gp,$x,$y,$size,$feed)=@_; 

  return if ($size<0); 
  my $holesize=$cp->{holesize};
  $gp->gcomment("Positioning for Hole"); 
  $gp->gmove('x',$x,'y',$y,'f',$feed); 
  
  my $z=0; 
  
  if ($holesize>$cp->{cuttersize}) 
  {
     $holesize-=$cp->{cuttersize}; # because we want to compensate for the size of the tool. 
     $gp->gmove('x',$x+$holesize/2,'y',$y); 
       
     my $passes=0; 
     while ($passes++ < $cp->{passes})
     {
        $z+=$cp->{passdepth}; # passdepth negative
        $gp->garcccw('x',$x-$holesize/2,'y',$y,'r',$holesize/2,'z',$z,'f',$feed); 
        $gp->garcccw ('x',$x+$holesize/2,'y',$y,'r',$holesize/2); 
     } 
     $gp->garcccw('x',$x-$holesize/2,'y',$y,'r',$holesize/2,'z',$z); # we always redo this as z depth was being faded in during this arc. 
  }
  else
  {
     $z=$cp->{passdepth}*$cp->{passes}; 
     $gp->gmove('z',$z,'f',$feed); 
  } 
  $gp->gmove('z', 0.1,$xs,$ys); 
  $gp->gcomment("Hole done");
}



sub cutwheel
{
#   my ($x,$y,$z,
#       $m,$np,$nw,
#       $gr,$cp,$dd,$dw,
#       $dp,$pf,$ad,$ar,$feed)=@_; 
    
    # 1 cut dededum.
    
    my ($cp ,   # wheel 
        $gp,    # graphics package, either generate graphics or gcode
        $x,$y,   # where to put the wheel 
        $feed
        )=@_;     

    my $t=0; # theta, angle of wheel;
    my $ti=0.5*360/$cp->{n}; # half tooth increment.   
    $ti*= $pi/180;      # in radians now. 

    
 
                        # In some situations particularly pinions tooth and gap angles are not the same. 
                        # define twf as factor extra for tooth, less than 1 for a wider gap 
    $tig=$ti*(2-$cp->{twf});  #  width of a gap in radians 
    $ti=$ti*$cp->{twf};       # this is now the width of a tooth. $tig+$ti is unchanged bu changes to $twf
    

    my ($xs,$ys)=($x,$y);
    
    $gp->ginit(); 
    $gp->gmove('z',0.1,'f',$feed); 

    $cp->cuthole($gp,$x,$y,$cp->{holesize},$feed);
    $cp->cuttrepan($gp,$xs,$ys);

    my $qtc=0.5*$cp->{mm}*$cp->{dw}*$pi/$cp->{n}; # quarter tooth circumference. 
    $x+=$cp->{mm}*$cp->{dw}*0.5+$qtc;
    $y+=$qtc; 
    $gp->gcomment("Move Away From Work"); 
    $gp->gmove('x',$x,'y',$y); 

    $gp->gcomment("Start Cutting"); 
    $gp->gmove('z',$cp->{passdepth}); 
    
    $x-= $qtc;
    $y-= $qtc; 
    $gp->gcompr('d',1,$gp->garccw('x',$x,'y',$y,'r',$qtc)); 

    $passes=0; 

    my $z=0; 
    while ($passes++ < $cp->{passes})
    {
    my $tcount=0;
    my $t=0; 
    $z+= $cp->{passdepth}; 
    while ($t/2.0/$pi<0.999 ) 
    {

    $gp->gcomment(sprintf("Tooth number %d pass %d ",++$tcount,$passes)); 
    
    $x-=$cp->{mm}*$cp->{dd}*cos($t);
    $y-=$cp->{mm}*$cp->{dd}*sin($t);
    # printf "G1 X$f Y$f Z$f F$ff\n", $x,$y,$z, $feed;    # radial stroke towards center of wheel
    
    $gp->gmove('x',$x,'y',$y,'z',$z,'f',$feed);
    
    $x+=$cp->{mm}*($cp->{dw}*0.5-$cp->{dd})*(cos($t+$tig)-cos($t));
    $y+=$cp->{mm}*($cp->{dw}*0.5-$cp->{dd})*(sin($t+$tig)-sin($t));
    # printf "G3 X$f Y$f R$f F$ff\n",$x,$y,$mm*($dw*0.5-$dd),$feed;      # bottom of tooth, flat bottom 
#    $gp->garccw('x',$x,'y',$y,'r',$cp->{mm}*($cp->{dw}-2*$cp->{dd})*$pi/$cp->{n}/4,'f',$feed); # semi circ bottom of tooth 
    $gp->garccw('x',$x,'y',$y,'r',$cp->{mm}*($cp->{dw}-2*$cp->{dd})*$tig/4,'f',$feed); # semi circ bottom of tooth 

    $t+=$tig;
    $x+=$cp->{mm}*$cp->{dd}*cos($t);         # radial stroke outwards
    $y+=$cp->{mm}*$cp->{dd}*sin($t); 
    $gp->gmove('x',$x,'y',$y); 

    # addendum 
    # add in addendum height.

    $x+=$cp->{mm}*$cp->{ad}*cos($t);
    $y+=$cp->{mm}*$cp->{ad}*sin($t); 
    # rotate to middle of tooth, 0.25 of full tooth+gap  width. 
    $x+=$cp->{mm}*($cp->{dw}*0.5+$cp->{ad})*(cos($t+0.5*$ti)-cos($t));
    $y+=$cp->{mm}*($cp->{dw}*0.5+$cp->{ad})*(sin($t+0.5*$ti)-sin($t));
    $gp->garcccw('x',$x,'y',$y,'r',$cp->{mm}*$cp->{ar},'f',$feed); 

    # back out the addendum height. 
    $x-=$cp->{mm}*$cp->{ad}*cos($t+0.5*$ti);
    $y-=$cp->{mm}*$cp->{ad}*sin($t+0.5*$ti); 
    # rotate a further half .25 tooth pitch 
    $x+=$cp->{mm}*($cp->{dw}*0.5)*(cos($t+$ti)-cos($t+0.5*$ti));
    $y+=$cp->{mm}*($cp->{dw}*0.5)*(sin($t+$ti)-sin($t+0.5*$ti));
    $gp->garcccw('x',$x,'y',$y,'r',$cp->{mm}*$cp->{ar},'f',$feed); 

    $t+=$ti; 

    }
    }

    $gp->gmove('z',0.1);
    $gp->gcomp0(); 
    $gp->gmove('x',$xs,'y',$ys); 

    $gp->gend(); 
}

sub dist
{
   my ($x1,$y1,$x2,$y2)=@_; 
   return sqrt(($x2-$x1)**2+($y2-$y1)**2);
}

sub smooth
{
    # given a circle radius b and a point on circumference l2
    # the plan is to smooth the join between the line l1/l2  (each of these are points) and the circle circumference by 
    # replacing a bit of the line l1/l2 with a circle of radius r. 
    # The following are supplied where x1 y1 is l1 point 1, x2, y2 point 2 or l2. 

    my ($x1,$y1,$x2,$y2,$b,$r)=@_;

    #print "smooth x1=$x1,y1=$y1,x2=$x2,y2=$y2,b=$b,r=$r\n";

    #($x1,$y1,$x2,$y2,$b,$r)=(-0.2 , 0.2 , -0.707 , 0.707 , 1.0 , 0.1); 
#    ($x1,$y1,$x2,$y2,$b,$r)=   (-0.279378067455943, 0.65017874253593, 0.702760341741972, 0.0831408677831007, 0.9,0.1); 

    my $ks; 
  
    # straight line l1 l2 has eqn y=mx+c 

    my $m=($y2-$y1)/($x2-$x1); 
    my $c=$y1-($y2-$y1)*$x1/($x2-$x1); 
    $ks=$r>0?1:-1;  
#   $ks=-$ks if ($y2>0);
    $r=abs($r);  
    $ks=-$ks if ($x1<$x2);

    
    # line paralell to this and distance r away is y=mx+j where j=c+k or j=c-k 

    my $k = abs($r) * sqrt( (($x2-$x1)**2+($y2-$y1)**2)/($x2-$x1)**2); 
    my $j=$c+$k*$ks;  # ks is the sign of k from above. 

    # we need to solve this with the circle x^2+y^2=(b+r)^2
    # substituting for y in here gives 
    # 
    # x^2+y^2=(b+r)^2
    # y=mx+j
    # x^2+m^2x^2+2mxj+j^2=(b+r)^2
    # (1+m^2) x^2 + 2mj x +j^2-(b+r)^2=0 
    # use quadratic equation formula to find x: 
    # x=(-b+- sqrt(b^2-4ac)/2a
    #
    # x=(-2mj +- sqrt(4m^2j^2-4(1+m^2)(j^2-(b+r)^2)))/2(1+m^2)
    # This is the center point of the arc. 

    # s is the distance between c and l2, we now have x and y coords for both c and l2. 
    # u^2=s^2-r^2 and is distance from l2 in direction of l1  for the new l2 point at start of smoothing arc. 
    # The end of the arc is oc scaled such that distance is b, the radius of the circle. 
    $r=-$r if (dist(0,0,$x1,$y1)<$b);
    
    my $cxa=(-2*$m*$j + sqrt(abs(4*$m**2*$j**2-4*(1+$m**2)*($j**2-($b+$r)**2))))/2/(1+$m**2);
    my $cxb=(-2*$m*$j - sqrt(abs(4*$m**2*$j**2-4*(1+$m**2)*($j**2-($b+$r)**2))))/2/(1+$m**2);   # This is the other root of the quadratic. 
                                                                                        # use the one closest to l2? 
    my $cya=$m*$cxa+$j; 
    my $cyb=$m*$cxb+$j; 

    ($cxa,$cya)=($cxb,$cyb) if (dist($x2,$y2,$cxb,$cyb)<dist($x2,$y2,$cxa,$cya)); # want nearest root to l2 . 

    my $s=dist($x2,$y2,$cxa,$cya);
    my $u=sqrt($s**2-$r**2);  
    my $sax=$x2+($x1-$x2)*$u/dist($x1,$y1,$x2,$y2);    #start of arc
    my $say=$y2+($y1-$y2)*$u/dist($x1,$y1,$x2,$y2); 
    my $eax=$cxa*$b/dist(0,0,$cxa,$cya); 
    my $eay=$cya*$b/dist(0,0,$cxa,$cya); 

    return ($x1,$y1,$sax,$say,$eax,$eay);              #ending on the circle

    # The code below is used fro graphing out test cases. 
    $gd=new gdcode("test.png",3.0,400,400); 
    $gd->gmove('z',0.1);
    $gd->gmove('x',$x1,'y',$y1); 
    $gd->gmove('z',-0.1); 
    $gd->gmove('x',$x2,'y',$y2); 

    $gd->gmove('z',0.1);         
    $gd->gmove('x',0,'y',$b); 
    $gd->gmove('z',-0.1);             
    $gd->garcccw('x',0,'y',-$b,'r',$b); 
    $gd->garcccw('x',0,'y',$b,'r',$b); 


    $gd->gmove('z',0.1); 
    $gd->gmove('x',$cxa+0.05,'y',$cya+0.05); 
    $gd->gmove('x',$cxa-0.05,'y',$cya-0.05,'z',-0.1); 
    $gd->gmove('x',$cxa+0.05,'y',$cya-0.05,'z',0.1); 
    $gd->gmove('x',$cxa-0.05,'y',$cya+0.05,'z',-0.1); 

#    $gd->gmove('x',$cxa,'y',$cya+$r,'z',0.1); 
#    $gd->garcccw('x',$cxa,'y',$cya-$r,'r',abs($r),'z',-0.1); 
#    $gd->garcccw('x',$cxa,'y',$cya+$r,'r',abs($r),'z',-0.1);

    
    $gd->gmove('x',$sax,'y',$say,'z',0.1); 
    $gd->garccw('x',$eax,'y',$eay,'r',abs($r),'z',-0.1);    


    $gd->gend(); 
    die; 
}
sub rsmooth
{
        my ($x1,$y1,$x2,$y2,$br,$r)=@_;  # params are point1, point 2, bossradius, radius
        #my @xy=swappair($x1,$y1,$x2,$y2); # point 1 and 2 given in wrong order for reverse smooth, so swap
        @xy=($x2,$y2,$x1,$y1); 
         print "rsmooth\n"; 
        @xy=smooth(@xy,$br,$r); # returns point start of arc, poin 1, start of arc end of arc
        @xy=@xy[4,5,2,3,0,1]; 
        # need in order given which is reverse, so give end of arc first
}


package cog; 

my  $pi=3.141592654;
my  $errorLimit=0.000001;
my $mm=1.0/24.8;  # 1mm is this inches;
my $inches="inches"; 
my $f="%9f  "; 
my $ff="%2.1f";   # for feed rate; 

sub newcogpair
{
    my ($this,$m,$np,$nw)=@_; 
    @_==4 or die "wrong number of paremeters";
    my $cogpair=bless {};

    my $w,$p; 

    $w=$cogpair->{wheel}=new wheel($m,$nw);
    $p=$cogpair->{pinion}=new wheel($m,$np);
    my $af=$cogpair->{af}=addendumFactor($np,$nw);

    $p->{pa}=$w->{pa}=$cogpair->{pa}=0.95*$af;  #practical addendum  factor
    $cogpair->{gr}=$np/$nw;#gear ratio
    $p->{cp}=$w->{cp}=$cogpair->{cp} = $m * $pi ;   # circular pitch
    $w->{dd}=$cogpair->{dd} = $m * $pi/2 ;
    $p->{dd}=$m*($af*0.95+0.4); # BSI rule for dd height for pinion. 
    $w->{dw}=$cogpair->{dw} = $m * $nw ;
    $p->{dw}=$cogpair->{dp} = $m * $np ;
#    $p->{ad}=$w->{ad}=$cogpair->{ad} = $m * 0.95 * $af ;
    $w->{ad}=$cogpair->{ad} = $m * 0.95 * $af ;
    $w->{ar}=$cogpair->{ar} = $m * 1.40 * $af ; 
    $w->{twf}=1.0;     # tooth width factor

    if ($p->{n}>=10)
    {
    # pinion profile A
    $p->{ar}=$m*0.525 ; 
    $p->{ad}=$m*0.525; 
    }
    elsif ($p->{n}==8 or $p->{n}==9)
    {
     # profile B. 
     $p->{ar}=$m*0.70 ; 
     $p->{ad}=$m*0.67; 
    }
    elsif ($p->{n}==6 or $p->{n}==7)
    {
      # profile C 
     $p->{ar}=$m*1.05; 
     $p->{ad}=$m*0.855; 
    }
    if ($p->{n}>=11) # set up special tooth width profile for pinion. 
    {
      $p->{twf}=1.25/1.57; 
    }
    else
    {
      $p->{twf}=1.05/1.57;       
    }

    $w->{mm}=$p->{mm}=$cogpair->{mm} = 1.0/24.8; 
    $cogpair->{nw}=$nw;
    $cogpair->{np}=$np; 

    return $cogpair; 
}

#compute(3.0,9,18);


sub addendumFactor
{   
    my ($np,$nw)=@_;
    my  $b  = 0.0 ;
    my  $t0 = 1.0 ;
    my  $t1 = 0.0 ;
    my  $r2 = 2 * $nw/$np ;
    while (abs($t1 - $t0) > $errorLimit)
        {   $t0 = $t1;
            $b = atan2(sin($t0), (1 + $r2 - cos($t0))) ;
            $t1 = $pi/$np + $r2 * $b ;  
        }
    return 0.25 * $np * (sin($t1)/sin($b) - $r2);
}





sub cutpinion_end
{
    my ($feed)=@_; 
    printf "G1 Z 0.1    F$ff\n", $feed; 
    print "G40\n"; 
}

sub cutpinion_init
{
             my $ti=0.5*360/$np; # half tooth increment.  
        $ti*= $pi/180;      # in radians now. 

        if ($np>=6 and $np<=10) # 6-10 teeth. 
        { 
          $toothwidthfactor=1.05/1.57; 
        } 
        elsif ($np>=11)
        {
          $toothwidthfactor=1.25/1.57; 
        }
       my $tit=$ti*$toothwidthfactor;  # tooth angular increment 
       my $tig=2*$ti-$tit;             # gap angular increment; 


        $x+=$mm*$dp/2+0.25*$mm*$tit*$dp;;
        $y+=0.25*$mm*$tit*$dp; 
        printf "G1 X$f Y$f\n" ,$x,$y; 
        printf "G1 Z$f\n",$z;

        $x-=0.25*$mm*$tit*$dp;
        $y-=0.25*$mm*$tit*$dp;
        printf "G1 Z$f F$ff\n",$z,$feed;
        printf "G42 D1 G2 X$f Y$f R$f F$ff\n",$x,$y,0.25*$mm*$tit*$dp,$feed;
        return ($x,$y); 
   

   
}


1;