(*:Name: FieldLines *)

(* :Title: FieldLines *)

(* :Author: Tom Wickham-Jones*)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 2.2 *)

(*:Summary:
	This package provides functions that calculate field lines.
*)

(* :History:
	Created summer 1993 by Tom Wickham-Jones.

	This package is described in the book
	Mathematica Graphics: Techniques and Applications.
	Tom Wickham-Jones, TELOS/Springer-Verlag 1994.

*)



BeginPackage[ "ExtendGraphics`FieldLines`",
				"Graphics`Arrow`"]
				
FieldLine::usage = 
	"FieldLine[ {x, ex, x0}, {y, ey, y0}, {t, t1}] will calculate
	the field line from the field {ex, ey}, starting at {x0,y0},
	of length t1.  x and y are the variables of the field and
	t is the variable of length down the trajectory."
	
FieldLine3D::usage = 
	"FieldLine3D[ {x,ex,x0}, {y,ey,y0}, {z,ez,z0}, {t, t1}] will calculate
	the field line from the field {ex,ey,ez}, starting at {x0,y0,z0},
	of length t1.  x, y and z are the variables of the field and
	t is the variable of length down the trajectory."
	
AddArrow::usage =
	"AddArrow[ line, dist] will travel down the line primitive
	and place arrows every dist for a maximum of 5 arrows.
	AddArrow[ line, dist, num] will generate num arrows."

Begin[ "`Private`"]


OnMessQ[ m_] := Head[m] =!= $Off

FieldLine[ {x_, ex_, x0_}, 
           {y_, ey_, y0_}, {t_, t1_}] :=
  Module[{exf, eyf, sol, t2, offM},
    {exf, eyf} = {ex, ey} /. {x -> x[t], y -> y[t]} ;
	offM = OnMessQ[ NDSolve::ndsz] ;
	If[ offM, Off[ NDSolve::ndsz]] ;
    sol = NDSolve[
            {x'[t] == exf, y'[t] == eyf,
             x[0] == x0, y[0] == y0
            }, {x[t], y[t]}, {t,0,t1}] ;
	If[ offM, On[ NDSolve::ndsz]] ;
    sol = {x[t], y[t]} /. First[ sol] ;
    If[ VectorQ[ sol /. t -> 0, NumberQ],
        t2 = Part[ sol, 1, 0, 1, 1, 2] ;
        sol =
            ParametricPlot[ Evaluate[ sol], {t,0,t2},
                DisplayFunction -> Identity] ;
  	  	First[ Cases[ sol, Line[_], Infinity]],
  	  	Line[ {{x0, y0}}]]
    ]


FieldLine[ {x_, ex_InterpolatingFunction, x0_}, 
           {y_, ey_InterpolatingFunction, y0_}, 
           {t_, t1_}] :=
  Module[{exf, eyf, x1, y1, sol, t2, oldM},
    exf[ x1_, y1_] := ex[ x1, y1] ;
    eyf[ x1_, y1_] := ey[ x1, y1] ;  
	offM = OnMessQ[ NDSolve::ndsz] ;
	If[ offM, Off[ NDSolve::ndsz]] ;
    sol =
      NDSolve[
        {x'[t] == exf[ x[t],y[t]], 
         y'[t] == eyf[ x[t],y[t]],
         x[0] == x0, y[0] == y0
        }, {x[t], y[t]}, {t,0,t1}] ;
	If[ offM, On[ NDSolve::ndsz]] ;
    sol = {x[t], y[t]} /. First[ sol] ;
    If[ VectorQ[ sol /. t -> 0, NumberQ],
        t2 = Part[ sol, 1, 0, 1, 1, 2] ;
        sol =
            ParametricPlot[ Evaluate[ sol], {t,0,t2},
                DisplayFunction -> Identity] ;
        First[ Cases[ sol, Line[_], Infinity]],
        Line[ {{x0, y0}}]]
    ]


(*

If this is to be used then the TriangularInterpolate
package should be loaded.  At present it is too slow.

FieldLine[ {x_, ex_TriangularInterpolating, x0_}, 
           {y_, ey_TriangularInterpolating, y0_}, 
           {t_, t1_}] :=
  Module[{exf, eyf, x1, y1, sol, t2},
    exf[ x1_, y1_] := ex[ x1, y1] ;
    eyf[ x1_, y1_] := ey[ x1, y1] ;  
    sol =
      NDSolve[
        {x'[t] == exf[ x[t],y[t]], 
         y'[t] == eyf[ x[t],y[t]],
         x[0] == x0, y[0] == y0
         }, {x[t], y[t]}, {t,0,t1}] ;
    sol = {x[t], y[t]} /. First[ sol] ;
    If[ VectorQ[ sol /. t -> 0, NumberQ],
        t2 = Part[ sol, 1, 0, 1, 1, 2] ;
        sol =
            ParametricPlot[ Evaluate[ sol], {t,0,t2},
                DisplayFunction -> Identity] ;
        First[ Cases[ sol, Line[_], Infinity]],
        Line[ {{x0, y0}}]]
    ]

*)

FieldLine3D[ {x_, ex_, x0_}, {y_, ey_, y0_}, 
             {z_, ez_, z0_}, {t_, t1_}] :=
  Module[{exf, eyf, ezf, sol, t2, oldM},
    {exf, eyf, ezf} = {ex, ey, ez} /. 
                   {x -> x[t], y -> y[t], z -> z[t]} ;
	offM = OnMessQ[ NDSolve::ndsz] ;
	If[ offM, Off[ NDSolve::ndsz]] ;
    sol =
      NDSolve[
        {
        x'[t] == exf, y'[t] == eyf, z'[t] == ezf,
        x[0] == x0, y[0] == y0, z[0] == z0
        }, {x[t], y[t], z[t]}, {t,0,t1}] ;
	If[ offM, On[ NDSolve::ndsz]] ;
    sol = {x[t], y[t], z[t]} /. First[ sol] ;
    If[ VectorQ[ sol /. t -> 0, NumberQ],
        t2 = Part[ sol, 1, 0, 1, 1, 2] ;
        sol =
            ParametricPlot3D[ Evaluate[ sol], {t,0,t2},
                DisplayFunction -> Identity] ;
        Cases[ sol, Line[_], Infinity],
        Line[ {{x0, y0, z0}}]]
    ]

AddArrow[ Line[ opts_], d_, num_:5] :=
  Module[{arrow = {}, n = 0, pts = Chop[ opts]},
    Fold[ 
      If[ First[ #1] >= d && n < num,
        n++;
        AppendTo[ arrow, 
                  Arrow[ Last[#1], #2, 
                         HeadScaling -> Absolute,
                         HeadCenter -> 0.5,
                         HeadLength -> 4]];
       {0, #2}, 
	   {First[#1] + 
        Sqrt[ Apply[ Plus, (Last[#1]-#2)^2]], #2}]&,
       {0, First[ pts]}, Rest[ pts]] ;
    arrow]				


End[]
				
				
			
EndPackage[]


				
				
				
				

