(**
		 DrawLine.m

(utility for FeynDraw)
**)

BeginPackage["DrawLine`","FeynArts`"]

DrawLine::usage = "
DrawLine[ from, to, center, particle, scale ] draws a propagator \n
line of particle. The line is straight for center=0 and part of \n
a circle for center={xcoord,ycoord}. scale is the total width of the \n
graph."

Orientation::usage = "
Orientation[ {x1,y1}, {x2,y2} ] gives the angle of this line against \n
the x-axis."

Begin["`Private`"]
                                   (** orientation of baseline **)

Orientation[p1_,p2_]:=Module[{xy1,xy2,phi},
  xy1=N[p1];
  xy2=N[p2];
  phi=0;
  If[ Abs[xy1[[1]]-xy2[[1]]]<10^-5,
          If[ xy2[[2]]-xy1[[2]] < 0, phi=-Pi/2, phi=Pi/2 ],
      phi=ArcTan[(xy2[[2]]-xy1[[2]])/(xy2[[1]]-xy1[[1]])];
      If[ Abs[xy2[[2]]-xy1[[2]]]<10^-5,phi=0 ];
      If[ xy2[[1]]<xy1[[1]],phi+=Pi ];
      If[ N[phi] > N[Pi] ,phi-= 2 Pi ]
    ];
  N[phi,10]          ]
				  (** constructing an arrow **)

Arrow[ X_, Y_, linedirection_, orientation_, arrowlength_ ] :=
  Module[ {alpha1,alpha2},
	 alpha1 = linedirection + Pi*(Sign[orientation]+1)/2 - Pi/6 ;
	 alpha2 = linedirection - Pi*(Sign[orientation]+1)/2 + Pi/6 ;
	 { Line[ {{ X+arrowlength*Cos[alpha1],
		    Y+arrowlength*Sin[alpha1] },
		  { X, Y },
		  { X+arrowlength*Cos[alpha2],
		    Y+arrowlength*Sin[alpha2] } } ] }
       ]
				  (** info about line form **)

GetLineSpec[ i_Integer ] := { straight,none,ToString[Global`c[i]] }

GetLineSpec[ part:type_[number__] ] := 
  If[ ValueQ[ LineSpec[ part ] ],
      Which[ Length[LineSpec[part]]==3,
	     LineSpec[part],
	     Length[LineSpec[part]]==2,
	     Join[ LineSpec[part],
	           { StringJoin@@ ToString/@ {type,number} } ],
	     Length[LineSpec[part]]==1,
	     Join[ LineSpec[part],
		   { none ,StringJoin@@ ToString/@ {type,number} } ]
	   ],
      {straight,none,StringJoin@@ ToString/@ {type,number} }
    ]

GetLineFlags[ { form_, direction_, symbol_ } ] :=
  Module[ {flaglist},
	flaglist = {symbol};
	Which[ direction===backward, flaglist=Prepend[ flaglist, -1 ],
	       direction===none    , flaglist=Prepend[ flaglist, 0 ],
	       direction===forward , flaglist=Prepend[ flaglist, 1 ] ];
	Which[ form===straight, flaglist=Join[ {0,0,0,0}, flaglist ],
	       form===wavy    , flaglist=Join[ {1,0,0,0}, flaglist ],
	       form===cycloid , flaglist=Join[ {1,1,0,0}, flaglist ],
	       form===dashed  , flaglist=Join[ {0,0,1,0}, flaglist ],
	       form===dotted  , flaglist=Join[ {0,0,0,1}, flaglist ] ];
	Return[ flaglist ]
       ]

DrawLine[ {Xa_,Ya_}, {Xb_,Yb_}, center_, particle_, 
          scaling_ ,thick_ ] :=
 Module[ {Xc,Yc,Xm,Ym,              (* center of line and circle *)
         lengthab,lengthac,         (*                           *)
	 height1,height2,radius,    (*                           *)
	 omegaab,omegaac,omegamc,   (* angles versus x-axis      *)
	 alpha,                     (* piechartangle             *)
	 numPi = N[ Pi, 10 ],       (*                           *)
	 drawpoints,                (* number of points          *)
	 amplitude,phaseamp,        (* real size of amplitudes   *)
	 wavenumber,                (* number of waves           *)
	 deltaX,                    (* x-increment               *)
	 deltaalpha,                (*                           *)
	 anglarray1,anglarray2,     (* arrays of angles and radii*)
	 radarray,XY1,XY2,          (*  -> array of points       *)
	 graphprim1,graphprim2,     (* graphics primitives       *)
	 graphprim3,                (*  (output)                 *)
	 flAmpl,flPhas,flDash,      (* flags                     *)
	 flDots,flArro,flTad,       (* flags                     *)
	 clockdir,                  (* clockwise/counterclockwise*)
         symb },                    (* symbol for propagator     *)
       Apply[ Set, #]& /@ 
	Transpose[{{flAmpl,flPhas,flDash,flDots,flArro,symb},
       		   GetLineFlags[ GetLineSpec[ particle ] ]    }];
(*
Print[" amplitude : ",flAmpl ];
Print[" phase     : ",flPhas ];
Print[" dashing   : ",flDash ];
Print[" dotted    : ",flDots ];
Print[" direction : ",flArro ];
Print[" symbol    : ",symb ];
*)
       lengthab = N[ Sqrt[ (Xa-Xb)^2+(Ya-Yb)^2 ], 10 ];
       If[ {Xa,Ya}==={Xb,Yb},
	   omegaab = 0. ; flTad=True,
           omegaab  = Orientation[ {Xa,Ya}, {Xb,Yb} ];flTad=False ];
(*
Print[" lengthab = ",lengthab];
Print[" omegaab  = ",180/numPi omegaab," deg"];
*)
       If[ center===0, 
	   If[ {Xa,Ya}==={Xb,Yb},
	       {Xc,Yc} = {Xa,Ya+2},
	       {Xc,Yc} = {N[(Xa+Xb)/2], Yc = N[(Ya+Yb)/2]} ],
	   {Xc,Yc} = {N[ center[[1]] ], N[ center[[2]] ]}  ];
       lengthac = N[ Sqrt[(Xa-Xc)^2+(Ya-Yc)^2], 10 ];
       omegaac  = Orientation[ {Xa,Ya}, {Xc,Yc} ];
       If[ omegaac > numPi ,omegaac -= 2 numPi  ];
       If[ omegaac <= - numPi ,omegaac += 2 numPi  ];
       If[ omegaac>-numPi/2 && omegaac<numPi/2 ,
           clockdir = Sign[ omegaac-omegaab ],
           clockdir = Sign[ omegaac-omegaab ]*
		      Sign[omegaac]*Sign[omegaab] ];
       If[ clockdir === 0 , clockdir = -1 ];
       omegamc  = If[ center===0, 
                      If[ flTad, omegaac, omegaab-(numPi/2)],
		      If[ flTad, 
			  omegaac,
                          If[omegaac>-numPi/2 && omegaac<numPi/2,
                             omegaab+ Sign[omegaac-omegaab] (numPi/2),
                             omegaab+ Sign[omegaac-omegaab] *
                                Sign[omegaac]*Sign[omegaab]*(numPi/2)] ]
		   ];
       If[ omegamc > numPi ,omegamc -= 2 numPi  ];
       If[ omegamc <= - numPi ,omegamc += 2 numPi  ];
(*
Print[" lengthac = ",lengthac];
Print[" omegaac  = ",180/numPi omegaac," deg"];
Print[" omegamc  = ",180/numPi omegamc," deg"];
*)
       If[ center===0 && !flTad ,
	   height1 = 0.;
	   radius  = 100.*scaling;
	   alpha   = N[ 2 ArcSin[ lengthab / 2 / radius ], 10 ];
	   height2 = N[ Abs[lengthab / 2 / Tan[ alpha/2 ]], 10 ],
	   height1 = N[ Abs[lengthac *
                             Sin[ArcCos[lengthab/2/lengthac]] ],10];
	   alpha   = If[ lengthab==0.,
			 2 numPi,
			 N[ 4 ArcTan[ 2 height1 / lengthab ], 10] ];
	   radius  = N[ lengthab / 2 / Sin[ alpha/2 ] , 10];
	   height2 = N[ radius Cos[ alpha/2 ], 10 ]
	 ];
       If[ flTad ,
	   {Xm,Ym}={(Xa+Xc)/2.,(Ya+Yc)/2.};
            radius=lengthac/2.;alpha= 2 numPi ,
           {Xm,Ym}={Xc,Yc}
		   +(height1+height2){-Cos[omegamc],-Sin[omegamc]}] ;
(*
Print[" alpha  = ",180/numPi alpha," deg"];
*)
       amplitude = flAmpl * 0.25 ;
       drawpoints = Floor[ If[ flAmpl == 0,
			   60 radius/scaling Abs[ alpha ],
			   100 radius/scaling Abs[alpha] (1+flPhas) ]];
       deltaX = lengthab / drawpoints;
       wavenumber = If[ flPhas == 1,
                    2 numPi (
		    Max[ 1, Floor[7.*Abs[alpha]*radius/scaling]]+.1),
		    numPi (
		    Max[ 1, Floor[10.*Abs[alpha]*radius/scaling]]+.5)];
       phaseamp = 0.8*flPhas * numPi/90./radius*scaling * Sign[alpha];
       (* now computing the points : *)
       deltaalpha = Array[ ( (#-1)/drawpoints * (alpha-phaseamp)/2.
			    )& , drawpoints+1 ];
       anglarray1 = Array[ ( omegamc - deltaalpha[[#]] - (phaseamp*
			     Sin[(#-1)/drawpoints*wavenumber])
			    )& , drawpoints+1 ];
       anglarray2 = Array[ ( omegamc + deltaalpha[[#]] + (phaseamp*
			     Sin[(#-1)/drawpoints*wavenumber])
			    )& , drawpoints+1 ];
       radarray = Array[ ( radius - flPhas * 1.5 * amplitude +
			   amplitude * (1+flPhas*0.5) * Cos[(#-1)*
			   wavenumber/drawpoints]
			  )& , drawpoints+1 ];
       XY1 = Array[ {Xm+radarray[[#]]*Cos[anglarray1[[#]]],
		     Ym+radarray[[#]]*Sin[anglarray1[[#]]]}& , 
		     drawpoints+1];
       XY2 = Array[ {Xm+radarray[[#]]*Cos[anglarray2[[#]]],
		     Ym+radarray[[#]]*Sin[anglarray2[[#]]]}& , 
		     drawpoints+1];
       (* now creating the output : *)
       graphprim1 = { Thickness[0.01*thick/scaling], Line[XY1] };
       graphprim2 = { Thickness[0.01*thick/scaling],Line[XY2] };
       Which[ flDash==1 ,
	      graphprim1 = Prepend[ graphprim1, Dashing[{.01,.01}] ];
	      graphprim2 = Prepend[ graphprim2, Dashing[{.01,.01}] ],
	      flDots==1 ,
	      graphprim1 = Prepend[ graphprim1, Dashing[{.002,.006}] ];
	      graphprim2 = Prepend[ graphprim2, Dashing[{.002,.006}] ]
	    ];
       If[ Head[symb]===Sub,
           graphprim3 = { Text[ FontForm[symb[[1]],Paint`Private`OutFont],
	                        {XY1[[1,1]] + 24./scaling Cos[omegamc], 
	                         XY1[[1,2]] + 24./scaling Sin[omegamc]} ],
		          Text[ FontForm[symb[[2]],Paint`Private`SmallOutFont],
			        {XY1[[1,1]] + 24./scaling (Cos[omegamc]+0.8), 
	                         XY1[[1,2]] + 24./scaling (Sin[omegamc]-0.5)}]
                        }, 
           graphprim3 = { Text[ FontForm[symb,Paint`Private`OutFont], 
                                {XY1[[1,1]] + 24./scaling Cos[omegamc],
			         XY1[[1,2]] + 24./scaling Sin[omegamc]}
                              ] } 
	 ];
       If[ flArro==0,
           Return[ {graphprim1,graphprim2,graphprim3} ],
	   Return[ {graphprim1,graphprim2,graphprim3,
                    Prepend[
		    Arrow[ N[ XY1[[1,1]] ], N[ XY1[[1,2]] ],
			   If[ center===0,
                               omegamc +  numPi/2 ,
			       omegamc - clockdir*numPi/2 ],
			   flArro ,
			   12./scaling ], 
                    Thickness[0.01*thick/scaling] ]                   } ]
(*                    Thickness[0.04/scaling] ]                   } ]*)
	]
]

End[]

EndPackage[]
