(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialiation Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)

BeginPackage["SegAnalysis`Two`"]
		
TwoA::usage =
  "TwoA[x, x1, t1, x2, t2] gives the expected fraction of A genotype at x, \
when two segregation distortion genes exist on the same gamete.  One trnsmits \
A genotype by t1 at the gene positition x1, Another transmits A genotype by \
t2 at the gene position x2."

TwoB::usage =
  "TwoB[x, x1, t1, x2, t2] gives the expected fraction of B genotype at x, \
when two segregation distortion genes exist on the same gamete.  One trnsmits \
A genotype by t1 at the gene positition x1, Another transmits A genotype by \
t2 at the gene position x2."

Begin["`Private`"]

(* === Warning messages ====  *)

SegregAnalysis::lmincrease = 
  "Warning: after `` decreases in the step size, NonlinearFit failed to find \
a new parameter estimate that will decrease the sum of squares."

(* ============   two    ================= *)
(* This uses Guess-Newton method for finding the least-
    squares solution for segregation distortion model by two genes on the 
      same gamete *)

mnlr[nptsA_, nptsB_, nresponseA_, nresponseB_, nresponseH_, x1_, 
   t1_, x2_, t2_, fp_, maxits_, wp_, pg_, ag_, mdls_, zt_] := 
  Module[{predA, predB, chilistA, chilistB, tmpchilistA, 
    tmpchilistB, tmpchilist, tmpnptsA, tmpnptsB, chisq, 
    derAx1, derBx1, derx1, derAt1, derBt1, dert1, derAx2, 
    derBx2, derx2, derAt2, derBt2, dert2, der, alpha, beta, 
    chisq0, newchi, tmpchi, its, accDelta, NAcc, precDelta, 
    NPrec, invaliddigits, accflag0, accflag1 = False, 
    accflag2 = False, precflag0, precflag1 = False, 
    precflag2 = False, ovflag = False, guessx1 = x1, 
    guesst1 = t1, guessx2 = x2, guesst2 = t2, oldx1, oldt1, 
    oldx2, oldt2, delta, deltax1, deltat1, deltax2, deltat2, 
    tmpx1, tmpt1, tmpx2, tmpt2, dits, dmaxits, newpredA, 
    newpredB, newchilistA, newchilistB, 
    accgoal = N[10^(-ag)], precgoal = N[10^(-pg)], 
    paramerr, errorDOF},
		(* calculation chisq *) 
		{{chilistA, chilistB}, chisq} = 
      calcchi[{nptsA, nptsB}, {nresponseA, nresponseB, nresponseH}, {guessx1, 
          guesst1, guessx2, guesst2}, wp];
    tmpchi = chisq; newchi = -1; its = 1;
		oldx1 = guessx1;oldt1 = guesst1;
		oldx2 = guessx2; oldt2 = guesst2;
		
		(* main loop of two *)	 
    ovflag = 
     Check[
				While[(
				(accDelta = Abs[chisq - newchi]; 
           NAcc = -Floor[Log[10, accDelta]] - 1; 
           precDelta = Abs[chisq - newchi]/(newchi + chisq)/2; 
           NPrec = -Floor[Log[10, precDelta]] - 1; 
           invaliddigits = TrueQ[
						NAcc =!= Infinity && 
              Accuracy[accDelta] < NAcc && 
              NPrec =!= Infinity && 
              Precision[precDelta] < Nprec]; 
           !invaliddigits ) && 
						
        ( accflag0 = ( 0 <= accDelta < accgoal) ;
		       precflag0 =( 0 <= precDelta < precgoal); 
           !((accflag0 && accflag1 && accflag2) || 
                (precflag0 && precflag1 && precflag2))  ) && 
        its <= maxits), 
					
			{accflag2, precflag2} = {accflag1, precflag1}; 
       {accflag1, precflag1} = {accflag0, precflag0}; 
        chisq = tmpchi; oldx1 = guessx1; oldt1 = guesst1; 
        oldx2 = guessx2; oldt2 = guesst2;
				(* find deltax1, deltat1, deltax2, 
            and deltat2 using LineaSolve for narmal equation *)	 
        tmpchilistA = Delete[chilistA, Position[nptsA, guessx1]]; 
        tmpnptsA = Delete[nptsA, Position[nptsA, guessx1]]; 
        tmpchilistB =Delete[chilistB, Position[nptsB, guessx1]]; 
        tmpnptsB =Delete[nptsB, Position[nptsB, guessx1]];
        tmpchilistA = Delete[tmpchilistA, Position[tmpnptsA, guessx2]]; 
        tmpnptsA = Delete[tmpnptsA, Position[tmpnptsA, guessx2]]; 
        tmpchilistB =Delete[tmpchilistB, Position[tmpnptsB, guessx2]]; 
        tmpnptsB =Delete[tmpnptsB, Position[tmpnptsB, guessx2]];
			 tmpchilist  = Join[tmpchilistA, tmpchilistB];
						(* calculate jacobian *)
					der={};
					delta={};
					If[!(MemberQ[fp, "x1"]),
						derAx1 =N[(Which[#1 <= guessx1, 1/2*
   (Sech[2*guessx1 - 2*#1]^2 - 
     2*guesst1*Sech[2*guessx1 - 2*#1]^2), 
  guessx1 < #1 && #1 < guessx2, 
  -((2*(-1 + Cosh[4*(guessx2 - #1)] - 
         2*Sinh[4*(guessx2 - #1)])*
       ((-5 + 18*guesst1 - 8*guesst2)*
          Cosh[4*(guessx1 - guessx2)] + 
         (-1 + 2*guesst1)*
          Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
         2*(1 + 2*guesst1 - 4*guesst2 + 
            Cosh[4*(guessx1 - #1)] + 
            2*guesst1*Cosh[4*(guessx1 - #1)] - 
            4*guesst2*Cosh[4*(guessx1 - #1)] + 
            Cosh[4*(guessx2 - #1)] + 
            2*guesst1*Cosh[4*(guessx2 - #1)] - 
            4*guesst2*Cosh[4*(guessx2 - #1)] - 
            2*Sinh[4*(guessx1 - guessx2)] + 
            8*guesst1*Sinh[4*(guessx1 - guessx2)] - 
            4*guesst2*Sinh[4*(guessx1 - guessx2)] + 
            2*Sinh[4*(guessx1 - #1)] - 
            4*guesst2*Sinh[4*(guessx1 - #1)] - 
            2*Sinh[4*(guessx2 - #1)] + 
            4*guesst2*Sinh[4*(guessx2 - #1)])))/
     (-7*Cosh[4*(guessx1 - guessx2)] + 
        Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
        2*(1 + Cosh[4*(guessx1 - #1)] + 
           Cosh[4*(guessx2 - #1)] - 
           4*Sinh[4*(guessx1 - guessx2)]))^2), 
  #1 >= guessx2, 0] & ) /@ tmpnptsA, wp]; 
			 derBx1 = 
         N[(Which[#1 <= guessx1, -(1/2) + guesst1 + 
   1/2*Tanh[2*(guessx1 - #1)]^2 - 
   guesst1*Tanh[2*(guessx1 - #1)]^2, 
  guessx1 < #1 && #1 < guessx2, 
  (2*(-1 + Cosh[4*(guessx2 - #1)] - 
       2*Sinh[4*(guessx2 - #1)])*
     ((-5 + 18*guesst1 - 8*guesst2)*
        Cosh[4*(guessx1 - guessx2)] + 
       (-1 + 2*guesst1)*
        Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
       2*(1 + 2*guesst1 - 4*guesst2 + 
          Cosh[4*(guessx1 - #1)] + 
          2*guesst1*Cosh[4*(guessx1 - #1)] - 
          4*guesst2*Cosh[4*(guessx1 - #1)] + 
          Cosh[4*(guessx2 - #1)] + 
          2*guesst1*Cosh[4*(guessx2 - #1)] - 
          4*guesst2*Cosh[4*(guessx2 - #1)] - 
          2*Sinh[4*(guessx1 - guessx2)] + 
          8*guesst1*Sinh[4*(guessx1 - guessx2)] - 
          4*guesst2*Sinh[4*(guessx1 - guessx2)] + 
          2*Sinh[4*(guessx1 - #1)] - 
          4*guesst2*Sinh[4*(guessx1 - #1)] - 
          2*Sinh[4*(guessx2 - #1)] + 
          4*guesst2*Sinh[4*(guessx2 - #1)])))/
   (-7*Cosh[4*(guessx1 - guessx2)] + 
      Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
      2*(1 + Cosh[4*(guessx1 - #1)] + 
         Cosh[4*(guessx2 - #1)] - 
         4*Sinh[4*(guessx1 - guessx2)]))^2, #1 >= guessx2, 0]  & ) /@ 
                  tmpnptsB, wp];
						derx1 = Join[derAx1, derBx1];
						der=Append[der,derx1]];
					
					If[!(MemberQ[fp, "t1"]),
						 derAt1 = N[(Which[#1 <= guessx1, 1/2*(1 - Tanh[2*(guessx1 - #1)]),
										 guessx1 < #1 && #1 < guessx2,
										(2*Cosh[2*(guessx1 - #1)]*
     (Cosh[2*(guessx1 - #1)] + Sinh[2*(guessx1 - #1)])*
     (-1 + Cosh[4*(guessx2 - #1)] - 
       2*Sinh[4*(guessx2 - #1)]))/
   (7*Cosh[4*(guessx1 - guessx2)] - 
     Cosh[4*(guessx1 + guessx2 - 2*#1)] - 
     2*(1 + Cosh[4*(guessx1 - #1)] + 
        Cosh[4*(guessx2 - #1)] - 
        4*Sinh[4*(guessx1 - guessx2)])), #1 >= guessx2, 0] & ) /@ tmpnptsA, 
                wp]; 
        derBt1 = 
         N[(Which[#1 <= guessx1, 1/2*(-1 + Tanh[2*(guessx1 - #1)]), 
  guessx1 < #1 && #1 < guessx2, 
  (2*Cosh[2*(guessx1 - #1)]*
     (Cosh[2*(guessx1 - #1)] + Sinh[2*(guessx1 - #1)])*
     (-1 + Cosh[4*(guessx2 - #1)] - 
       2*Sinh[4*(guessx2 - #1)]))/
   (-7*Cosh[4*(guessx1 - guessx2)] + 
     Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
     2*(1 + Cosh[4*(guessx1 - #1)] + 
        Cosh[4*(guessx2 - #1)] - 
        4*Sinh[4*(guessx1 - guessx2)])), #1 >= guessx2, 0] & ) /@ tmpnptsB, 
                wp];
						dert1 = Join[derAt1, derBt1];
						der=Append[der,dert1]];
					
					If[!(MemberQ[fp, "x2"]),
						derAx2 = N[(Which[#1 <= guessx1, 0, guessx1 < #1 && #1 < guessx2, 
  -((2*(-1 + Cosh[4*(guessx1 - #1)] + 
         2*Sinh[4*(guessx1 - #1)])*
       ((5 + 8*guesst1 - 18*guesst2)*
          Cosh[4*(guessx1 - guessx2)] + 
         (1 - 2*guesst2)*
          Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
         2*(-1 + 4*guesst1 - 2*guesst2 - 
            Cosh[4*(guessx1 - #1)] + 
            4*guesst1*Cosh[4*(guessx1 - #1)] - 
            2*guesst2*Cosh[4*(guessx1 - #1)] - 
            Cosh[4*(guessx2 - #1)] + 
            4*guesst1*Cosh[4*(guessx2 - #1)] - 
            2*guesst2*Cosh[4*(guessx2 - #1)] + 
            2*Sinh[4*(guessx1 - guessx2)] + 
            4*guesst1*Sinh[4*(guessx1 - guessx2)] - 
            8*guesst2*Sinh[4*(guessx1 - guessx2)] - 
            2*Sinh[4*(guessx1 - #1)] + 
            4*guesst1*Sinh[4*(guessx1 - #1)] + 
            2*Sinh[4*(guessx2 - #1)] - 
            4*guesst1*Sinh[4*(guessx2 - #1)])))/
     (-7*Cosh[4*(guessx1 - guessx2)] + 
        Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
        2*(1 + Cosh[4*(guessx1 - #1)] + 
           Cosh[4*(guessx2 - #1)] - 
           4*Sinh[4*(guessx1 - guessx2)]))^2), 
  #1 >= guessx2, 1/2*(-1 + 2*guesst2)*
   Sech[2*guessx2 - 2*#1]^2] & ) /@ tmpnptsA, wp]; 
			 derBx2 = 
         N[(Which[#1 <= guessx1, 0, guessx1 < #1 && #1 < guessx2, 
  (2*(-1 + Cosh[4*(guessx1 - #1)] + 
       2*Sinh[4*(guessx1 - #1)])*
     ((5 + 8*guesst1 - 18*guesst2)*
        Cosh[4*(guessx1 - guessx2)] + 
       (1 - 2*guesst2)*
        Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
       2*(-1 + 4*guesst1 - 2*guesst2 - 
          Cosh[4*(guessx1 - #1)] + 
          4*guesst1*Cosh[4*(guessx1 - #1)] - 
          2*guesst2*Cosh[4*(guessx1 - #1)] - 
          Cosh[4*(guessx2 - #1)] + 
          4*guesst1*Cosh[4*(guessx2 - #1)] - 
          2*guesst2*Cosh[4*(guessx2 - #1)] + 
          2*Sinh[4*(guessx1 - guessx2)] + 
          4*guesst1*Sinh[4*(guessx1 - guessx2)] - 
          8*guesst2*Sinh[4*(guessx1 - guessx2)] - 
          2*Sinh[4*(guessx1 - #1)] + 
          4*guesst1*Sinh[4*(guessx1 - #1)] + 
          2*Sinh[4*(guessx2 - #1)] - 
          4*guesst1*Sinh[4*(guessx2 - #1)])))/
   (-7*Cosh[4*(guessx1 - guessx2)] + 
      Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
      2*(1 + Cosh[4*(guessx1 - #1)] + 
         Cosh[4*(guessx2 - #1)] - 
         4*Sinh[4*(guessx1 - guessx2)]))^2, #1 >= guessx2, 
  -(1/2)*(-1 + 2*guesst2)*Sech[2*guessx2 - 2*#1]^2] & ) /@ tmpnptsB, wp];
						derx2 = Join[derAx2, derBx2];
						der=Append[der,derx2]];
					
					If[!(MemberQ[fp, "t2"]),
						derAt2 =N[(Which[#1 <= guessx1, 0, guessx1 < #1 && #1 < guessx2, 
  (Cosh[2*(guessx2 - #1)]*
     (-1 + Cosh[4*(guessx1 - #1)] + 
       2*Sinh[4*(guessx1 - #1)])*
     (-2*Cosh[2*(guessx2 - #1)] + 2*Sinh[2*(guessx2 - #1)])
)/(-7*Cosh[4*(guessx1 - guessx2)] + 
     Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
     2*(1 + Cosh[4*(guessx1 - #1)] + 
        Cosh[4*(guessx2 - #1)] - 
        4*Sinh[4*(guessx1 - guessx2)])), #1 >= guessx2, 
  1/2*(1 + Tanh[2*(guessx2 - #1)])] & ) /@ tmpnptsA, wp]; 
        derBt2 = 
         N[(Which[#1 <= guessx1, 0, guessx1 < #1 && #1 < guessx2, 
  (4*Cosh[2*(guessx2 - #1)]*Sinh[2*(guessx1 - #1)]*
     (2*Cosh[2*(guessx1 - #1)] + Sinh[2*(guessx1 - #1)])*
     (Cosh[2*(guessx2 - #1)] - Sinh[2*(guessx2 - #1)]))/
   (-7*Cosh[4*(guessx1 - guessx2)] + 
     Cosh[4*(guessx1 + guessx2 - 2*#1)] + 
     2*(1 + Cosh[4*(guessx1 - #1)] + 
        Cosh[4*(guessx2 - #1)] - 
        4*Sinh[4*(guessx1 - guessx2)])), #1 >= guessx2, 
  1/2*(-1 + Tanh[2*(-guessx2 + #1)])] & ) /@ tmpnptsB, wp];
						dert2 = Join[derAt2, derBt2];
						der=Append[der, dert2]]; 

			(* normal equation: alpha . deltaguess === beta   *)
			 {alpha, beta} = 
            N[{der . Transpose[der], Map[tmpchilist . #1 &, der]}, wp];
					
					delta = N[LinearSolve[alpha, beta], wp];
					If[!(MemberQ[fp, "x1"]),deltax1=First[delta];delta=Rest[delta], 
            deltax1=0];
					If[!(MemberQ[fp, "t1"]),deltat1=First[delta]; delta=Rest[delta],  
            deltat1=0];
					If[!(MemberQ[fp, "x2"]),deltax2=First[delta];delta=Rest[delta], 
            deltax2=0];
					If[!(MemberQ[fp, "t2"]),deltat2=First[delta],  deltat2=0];
					
				 {tmpx1, tmpt1, tmpx2, tmpt2} = 
           
            N[{oldx1, oldt1, oldx2, oldt2} +  {deltax1, deltat1, deltax2, 
                  deltat2}, wp]; 
					{{newchilistA, newchilistB}, newchi} = 
            calcchi[{nptsA, nptsB}, {nresponseA, nresponseB, nresponseH}, {
                tmpx1, tmpt1, tmpx2, tmpt2}, wp];
			(* find newguess for smaller chisq  *) 				
			 dits = 1; dmaxits = 15; 
        While[newchi - chisq >= 0 && dits <= dmaxits, 
         {tmpx1, tmpt1, tmpx2, tmpt2} = 
           
              N[{oldx1, oldt1, oldx2, 
                    oldt2} + (0.5 ^ dits) {deltax1, deltat1, deltax2, 
                      deltat2}, wp]; 
				{{newchilistA, newchilistB}, newchi} = 
              calcchi[{nptsA, nptsB}, {nresponseA, nresponseB, nresponseH}, {
                  tmpx1, tmpt1, tmpx2, tmpt2}, wp];
				 dits++]; (* end While[newchi     *)
        If[dits > dmaxits, Message[SegregAnalysis::lmincrease, dmaxits]]; 
        tmpchi = newchi; 
        {guessx1, guesst1, guessx2, guesst2} = {tmpx1, tmpt1, tmpx2, tmpt2}; 
        chilistA = newchilistA; chilistB = newchilistB; 
        printprogress[its, newchi, 
         {guessx1, guesst1, guessx2, guesst2}];
		 its++], (* end While *)
      True, General::ovfl]; (* end Check *)
			(* Adjust the accuracy and precision of chisq *)
    If[Accuracy[newchi] > NAcc, 
     newchi = SetAccuracy[newchi, NAcc]]; 
    If[Precision[newchi] > NPrec, 
     newchi = SetPrecision[newchi, NPrec]]; 
    errorDOF = 
     Length[chilistA] + Length[chilistB] + Length[nresponseH] -4+Length[fp];
		paramerr = Sqrt[(newchi/errorDOF) DiagonalElements[Inverse[alpha]]]; 
    {newchi, {guessx1, guesst1, guessx2, guesst2}, paramerr, errorDOF, fp}
	]

(*  *)
TwoA[x_, x1_, t1_, x2_, t2_] :=
  Which[x <= x1, 1/2*(t1 + 1/2*Tanh[2*(-x + x1)] - 
     t1*Tanh[2*(-x + x1)]), x1 < x && x < x2, 
  (-4*t1*(Cosh[4*x1] + Sinh[4*x1])*
      (Cosh[4*x] + Cosh[4*x1] + Sinh[4*x] + Sinh[4*x1])*
      (-Cosh[8*x2] + 3*(Cosh[8*x] + Sinh[8*x]) - 
        Sinh[8*x2] - 2*(Cosh[4*(x + x2)] + 
           Sinh[4*(x + x2)])) + 
     (Cosh[4*x] - Cosh[4*x1] + Sinh[4*x] - Sinh[4*x1])*
      ((-1 + 4*t2)*(Cosh[12*x] + Sinh[12*x]) + 
        3*(1 + 4*t2)*(Cosh[4*(2*x + x1)] + 
           Sinh[4*(2*x + x1)]) + 
        (Cosh[4*x] + Cosh[4*x1] + Sinh[4*x] + Sinh[4*x1])*
         (Cosh[8*x2] + Sinh[8*x2]) + 
        4*(t2*(Cosh[4*x] + Sinh[4*x]) + 
           (-1 + 3*t2)*(Cosh[4*x1] + Sinh[4*x1]))*
         (Cosh[4*(x + x2)] + Sinh[4*(x + x2)])))/
   (4*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
       Sinh[8*x] - 3*(Cosh[4*(x + x1)] + 
          Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
       Sinh[4*(x1 + x2)])*
     (Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
       Sinh[8*x] + 5*(Cosh[4*(x + x1)] + 
          Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
       Sinh[4*(x1 + x2)])), x >= x2, 
  1/2*(t2 + 1/2*Tanh[2*(x - x2)] - t2*Tanh[2*(x - x2)])]
		
(*  *)
TwoB[x_, x1_, t1_, x2_, t2_] :=
  Which[x <= x1, (1 - t1 - Tanh[2*(-x + x1)]/2 + 
     t1*Tanh[2*(-x + x1)])/2, x1 < x && x < x2, 
  (4*t1*(Cosh[4*x1] + Sinh[4*x1])*
      (Cosh[4*x] + Cosh[4*x1] + Sinh[4*x] + Sinh[4*x1])*
      (-Cosh[8*x2] + 3*(Cosh[8*x] + Sinh[8*x]) - 
        Sinh[8*x2] - 2*(Cosh[4*(x + x2)] + 
           Sinh[4*(x + x2)])) - 
     (Cosh[4*x] + Sinh[4*x] + 3*(Cosh[4*x1] + Sinh[4*x1]))*
      (4*t2*(Cosh[4*x] + Sinh[4*x])*
         (Cosh[4*x] - Cosh[4*x1] + Sinh[4*x] - 
           Sinh[4*x1])*(Cosh[4*x] + Cosh[4*x2] + 
           Sinh[4*x] + Sinh[4*x2]) + 
        (-Cosh[4*x2] - 3*(Cosh[4*x] + Sinh[4*x]) - 
           Sinh[4*x2])*(Cosh[8*x] + Cosh[4*(x + x2)] + 
           Cosh[4*(x1 + x2)] + Sinh[8*x] - 
           3*(Cosh[4*(x + x1)] + Sinh[4*(x + x1)]) + 
           Sinh[4*(x + x2)] + Sinh[4*(x1 + x2)])))/
   (4*(Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
       Sinh[8*x] - 3*(Cosh[4*(x + x1)] + 
          Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
       Sinh[4*(x1 + x2)])*
     (Cosh[8*x] + Cosh[4*(x + x2)] + Cosh[4*(x1 + x2)] + 
       Sinh[8*x] + 5*(Cosh[4*(x + x1)] + 
          Sinh[4*(x + x1)]) + Sinh[4*(x + x2)] + 
       Sinh[4*(x1 + x2)])), x >= x2, 
  (1 - t2 - Tanh[2*(x - x2)]/2 + t2*Tanh[2*(x - x2)])/2]

(*   calculation of chi-list and chi-squares *)
calcchi[pts_List, response_List, guess_List, wp_] :=
	Block[{x1,t1, x2, t2, ptsA, ptsB, ptsH, responseA, responseB, responseH, 
      predA, predB, predH, chilistA, chilistB, chilistH, chisq},
		{responseA, responseB, responseH} = response;
		{x1, t1, x2, t2} = N[guess,wp];
		{ptsA, ptsB} = pts;
		predA = N[TwoA[#1, x1, t1, x2, t2] & /@ ptsA, wp];
		predB = N[TwoB[#1, x1, t1, x2, t2] & /@ ptsB, wp];
		chilistA = responseA - predA;
		chilistB = responseB - predB;
		chilistH = N[ (#1 - 0.5 )& /@ responseH, wp];
		chisq =N[ 
        Plus @@ (
            Join[Map[#1^2 &, chilistA], Map[#1^2 &, chilistB],  
              Map[#1^2 &, chilistH]]), wp];
		{{chilistA, chilistB}, chisq}
		]

(*    *)
printprogress[iteration_, chisq_, params_] := 
     Print[
    StringForm[
      "Iteration:`1` ChiSquared:`2`  {Position, Transmittance}:`3`", 
											NumberForm[iteration, 6], NumberForm[chisq, 6], 
      NumberForm[params, 6]]]

(*    *)
DiagonalElements[matrix_] := 
  Module[{k, size = Length[matrix]}, 
   Table[matrix[[k,k]], {k, 1, size}]
		]

(* =========================================== *)
End[]
EndPackage[]