(***********************************************************************
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`OneOne`"]

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

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

OneOneH::usage =
  "OneOneH[x, x1, t1, x2, t2] gives the expected fraction of H genotype at x, \
when two segregation distortion genes exit on the different 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."

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

mnlr[nptsA_, nptsB_, nptsH_, nresponseA_, nresponseB_, nresponseH_, x1_, t1_, 
    x2_, t2_, fp_, maxits_, 
   wp_, pg_, ag_, mdls_, zt_] := 
  Module[{predA, predB, predH, chilistA, chilistB, 
    chilistH, tmpchilistA, tmpchilistB, tmpchilistH, 
    tmpchilist, tmpnptsA, tmpnptsB, tmpnptsH, chisq, 
    derAx1, derBx1, derHx1, derx1, derAt1, derBt1, derHt1, 
    dert1, derAx2, derBx2, derHx2, derx2, derAt2, derBt2, 
    derHt2, 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, newpredH, newchilistA, newchilistB, 
    newchilistH, accgoal = N[10^(-ag)], 
    precgoal = N[10^(-pg)], paramerr, errorDOF}, 
 (* calculating chisq *)
		{{chilistA, chilistB, chilistH}, chisq} = 
      calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, nresponseH}, {
          guessx1, guesst1, guessx2, guesst2}, wp];
	 tmpchi = chisq; newchi = -1;its = 1;
		oldx1 = guessx1; oldt1 = guesst1;
		oldx2 = guessx2; oldt2 = guesst2; 
    
		(* main loop for oneone *)
		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]]; 
        tmpchilistH =Delete[chilistH, Position[nptsH, guessx1]]; 
        tmpnptsH =Delete[nptsH, Position[nptsH, 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]]; 
        tmpchilistH = Delete[tmpchilistH, Position[tmpnptsH, guessx2]]; 
        tmpnptsH = Delete[tmpnptsH, Position[tmpnptsH, guessx2]]; 
        tmpchilist=Join[tmpchilistA, tmpchilistB, tmpchilistH];
					(* calculate jacobian *)
					der={};
					delta={};
						If[!(MemberQ[fp, "x1"]),
						derAx1=N[(Which[#1 < guessx1, 
              guesst2*Sech[2*Abs[#1 - guessx1]]^2 - 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2 + 
               1/2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] - 
               guesst1*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] - 
               guesst2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] + 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]], 
              #1 > guessx1, 
              -guesst2*Sech[2*Abs[#1 - guessx1]]^2 + 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2 - 
               1/2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] + 
               guesst1*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] + 
               guesst2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] - 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]]] & ) /@            tmpnptsA, wp];
						derBx1 = N[(Which[#1 < guessx1, 
              -Sech[2*Abs[#1 - guessx1]]^2 + 
               2*guesst1*
                Sech[2*Abs[#1 - guessx1]]^2 + 
               guesst2*Sech[2*Abs[#1 - guessx1]]^2 - 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2 + 
               1/2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] - 
               guesst1*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] - 
               guesst2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] + 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]], 
              #1 > guessx1, 
              Sech[2*Abs[#1 - guessx1]]^2 - 
               2*guesst1*Sech[2*Abs[#1 - guessx1]]^2 - 
               guesst2*Sech[2*Abs[#1 - guessx1]]^2 + 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2 - 
               1/2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] + 
               guesst1*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] + 
               guesst2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] - 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]]] & ) /@ tmpnptsB, wp];
						derHx1 = N[(Which[#1 < guessx1, 
              Sech[2*Abs[#1 - guessx1]]^2 - 
               2*guesst1*
                Sech[2*Abs[#1 - guessx1]]^2 - 
               2*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2 + 
               4*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2 - 
               Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] + 
               2*guesst1*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] + 
               2*guesst2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] - 
               4*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]], 
              #1 > guessx1, 
              -Sech[2*Abs[#1 - guessx1]]^2 + 
               2*guesst1*Sech[2*Abs[#1 - guessx1]]^2 + 
               2*guesst2*Sech[2*Abs[#1 - guessx1]]^2 - 
               4*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2 + 
               Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] - 
               2*guesst1*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] - 
               2*guesst2*Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]] + 
               4*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx1]]^2*
                Tanh[2*Abs[#1 - guessx2]]] & ) /@ tmpnptsH, wp];
						derx1 = Join[derAx1, derBx1, derHx1];
						der=Append[der,derx1]];
					If[!(MemberQ[fp, "t1"]),
						derAt1 =N[(guesst2 - guesst2*
               Tanh[2*Abs[#1 - guessx1]] + 
              1/2*Tanh[2*Abs[#1 - guessx2]] - 
              guesst2*Tanh[2*Abs[#1 - guessx2]] - 
              1/2*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] + 
              guesst2*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] & ) /@ tmpnptsA, wp];
						derBt1 = N[(-1 + guesst2 + Tanh[2*Abs[#1 - guessx1]] - 
              guesst2*Tanh[2*Abs[#1 - guessx1]] + 
              1/2*Tanh[2*Abs[#1 - guessx2]] - 
              guesst2*Tanh[2*Abs[#1 - guessx2]] - 
              1/2*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] + 
              guesst2*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] & ) /@ tmpnptsB, wp];
						derHt1 =N[(1 - 2*guesst2 - 
              Tanh[2*Abs[#1 - guessx1]] + 
              2*guesst2*Tanh[2*Abs[#1 - guessx1]] - 
              Tanh[2*Abs[#1 - guessx2]] + 
              2*guesst2*Tanh[2*Abs[#1 - guessx2]] + 
              Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] - 
              2*guesst2*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] & ) /@ tmpnptsH, wp];
						dert1 = Join[derAt1, derBt1, derHt1];
						der=Append[der,dert1]];
					If[!(MemberQ[fp, "x2"]),
						derAx2 =N[(Which[#1 < guessx2, 
              guesst1*Sech[2*Abs[#1 - guessx2]]^2 - 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2 + 
               1/2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] - 
               guesst1*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] - 
               guesst2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] + 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]], 
              #1 > guessx2, 
              -guesst1*Sech[2*Abs[#1 - guessx2]]^2 + 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2 - 
               1/2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] + 
               guesst1*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] + 
               guesst2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] - 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]]] & ) /@ tmpnptsA, wp];
						derBx2 =N[(Which[#1 < guessx2, 
              -Sech[2*Abs[#1 - guessx2]]^2 + 
               guesst1*Sech[2*Abs[#1 - guessx2]]^2 + 
               2*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2 - 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2 + 
               1/2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] - 
               guesst1*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] - 
               guesst2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] + 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]], 
              #1 > guessx2, 
              Sech[2*Abs[#1 - guessx2]]^2 - 
               guesst1*Sech[2*Abs[#1 - guessx2]]^2 - 
               2*guesst2*Sech[2*Abs[#1 - guessx2]]^2 + 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2 - 
               1/2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] + 
               guesst1*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] + 
               guesst2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] - 
               2*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]]] & ) /@ tmpnptsB, wp];
						derHx2 = N[(Which[#1 < guessx2, 
              Sech[2*Abs[#1 - guessx2]]^2 - 
               2*guesst1*
                Sech[2*Abs[#1 - guessx2]]^2 - 
               2*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2 + 
               4*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2 - 
               Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] + 
               2*guesst1*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] + 
               2*guesst2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] - 
               4*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]], 
              #1 > guessx2, 
              -Sech[2*Abs[#1 - guessx2]]^2 + 
               2*guesst1*Sech[2*Abs[#1 - guessx2]]^2 + 
               2*guesst2*Sech[2*Abs[#1 - guessx2]]^2 - 
               4*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2 + 
               Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] - 
               2*guesst1*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] - 
               2*guesst2*Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]] + 
               4*guesst1*guesst2*
                Sech[2*Abs[#1 - guessx2]]^2*
                Tanh[2*Abs[#1 - guessx1]]] & ) /@ tmpnptsH, wp];
						derx2 = Join[derAx2, derBx2, derHx2];
						der=Append[der,derx2]];
					If[!(MemberQ[fp, "t2"]),
						derAt2 =N[(guesst1 + 1/2*Tanh[2*Abs[#1 - guessx1]] - 
              guesst1*Tanh[2*Abs[#1 - guessx1]] - 
              guesst1*Tanh[2*Abs[#1 - guessx2]] - 
              1/2*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] + 
              guesst1*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] & ) /@ tmpnptsA, wp];
						derBt2 =N[(-1 + guesst1 + 
              1/2*Tanh[2*Abs[#1 - guessx1]] - 
              guesst1*Tanh[2*Abs[#1 - guessx1]] + 
              Tanh[2*Abs[#1 - guessx2]] - 
              guesst1*Tanh[2*Abs[#1 - guessx2]] - 
              1/2*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] + 
              guesst1*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] & ) /@ tmpnptsB, wp];
						derHt2 =N[(1 - 2*guesst1 - 
              Tanh[2*Abs[#1 - guessx1]] + 
              2*guesst1*Tanh[2*Abs[#1 - guessx1]] - 
              Tanh[2*Abs[#1 - guessx2]] + 
              2*guesst1*Tanh[2*Abs[#1 - guessx2]] + 
              Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] - 
              2*guesst1*Tanh[2*Abs[#1 - guessx1]]*
               Tanh[2*Abs[#1 - guessx2]] & ) /@ tmpnptsH, wp];
						dert2 = Join[derAt2, derBt2, derHt2];
						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];
			(* calculation chi-squares *)
			{{newchilistA, newchilistB, newchilistH}, newchi} = 
            calcchi[{nptsA, nptsB, nptsH}, {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, newchilistH}, newchi} = 
              calcchi[{nptsA, nptsB, nptsH}, {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; 
        chilistH = newchilistH; 
        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[chilistH]- 4+Length[fp];
		paramerr = Sqrt[(newchi/errorDOF) DiagonalElements[Inverse[alpha]]]; 
    
 {newchi, {guessx1, guesst1, guessx2, guesst2}, paramerr, errorDOF, fp}
]

(*  *)
OneOneA[x_, x1_, t1_, x2_, t2_] :=
  t1*t2 + 1/2*t2*Tanh[2*Abs[x - x1]] - t1*t2*Tanh[2*Abs[x - x1]] + 
    1/2*t1*Tanh[2*Abs[x - x2]] - t1*t2*Tanh[2*Abs[x - x2]] + 
    1/4*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]] - 
    1/2*t1*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]] - 
    1/2*t2*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]] + 
    t1*t2*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]]

(*   *)
OneOneB[x_, x1_, t1_, x2_, t2_] :=
  1 - t1 - t2 + t1*t2 - 1/2*Tanh[2*Abs[x - x1]] + t1*Tanh[2*Abs[x - x1]] + 
    1/2*t2*Tanh[2*Abs[x - x1]] - t1*t2*Tanh[2*Abs[x - x1]] - 
    1/2*Tanh[2*Abs[x - x2]] + 1/2*t1*Tanh[2*Abs[x - x2]] + 
    t2*Tanh[2*Abs[x - x2]] - t1*t2*Tanh[2*Abs[x - x2]] + 
    1/4*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]] - 
    1/2*t1*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]] - 
    1/2*t2*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]] + 
    t1*t2*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]]

(*  *)
OneOneH[x_, x1_, t1_, x2_, t2_] :=
  t1 + t2 - 2*t1*t2 + 1/2*Tanh[2*Abs[x - x1]] - t1*Tanh[2*Abs[x - x1]] - 
    t2*Tanh[2*Abs[x - x1]] + 2*t1*t2*Tanh[2*Abs[x - x1]] + 
    1/2*Tanh[2*Abs[x - x2]] - t1*Tanh[2*Abs[x - x2]] - 
    t2*Tanh[2*Abs[x - x2]] + 2*t1*t2*Tanh[2*Abs[x - x2]] - 
    1/2*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]] + 
    t1*Tanh[2*Abs[x - x1]]* Tanh[2*Abs[x - x2]] + 
    t2*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]] - 
    2*t1*t2*Tanh[2*Abs[x - x1]]*Tanh[2*Abs[x - x2]]

(*   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, ptsH} = pts;
		predA = N[OneOneA[#1, x1, t1, x2, t2] & /@ ptsA, wp];
		predB = N[OneOneB[#1, x1, t1, x2, t2] & /@ ptsB, wp];
		predH = N[OneOneH[#1, x1, t1, x2, t2] & /@ ptsH, wp];
		chilistA = responseA - predA;
		chilistB = responseB - predB;
		chilistH = responseH - predH;
		chisq =N[ 
        Plus @@ (
            Join[Map[#1^2 &, chilistA], Map[#1^2 &, chilistB], 
              Map[#1^2&, chilistH]]), wp];
			{{chilistA, chilistB, chilistH}, 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[]