# Analytically finding critical values in a hysteresis curve made with ListPlot

Let me start with a super basic explanation of what I'm trying to do, and I will follow with a more detailed discussion afterwards.

From data in the form shown below (a Table[] of (x,y) pairs), I would like to

1. Extract the coordinates of the left-most and right-most points. I know how to find the minimum and maximum x-coordinates but not the corresponding y-coordinates.
2. Extract the x-value for which a sharp transition occurs between two saturated y-values.

I'm trying to figure out if there's a way to find things like saturation and coercivity in a hysteresis curve from real data using Mathematica. For those not familiar, hysteresis curves generally looks something like this:

Pardon the overlay of 3 different plots but this is particularly helpful as it shows the sometimes varied nature of such curves.

I have a variable which is a Table of (x,y) pairs of my data,

pairs =
Table[{xvalues[[x]], yvalues[[x]]}, {x, 1,
Length[xvalues]}]

From this data I would like to extract the saturation and coercivity (labeled M_s and H_c respectively in the image above). The curves are not in general not centered about H=0 and M=0.

A common definition for the coercivity is to take the x value for which the y value is half way between its min and max, and this is what I was planning on using. A possible issue is that most of my curves don't have a physical data point corresponding to this because the transition is very sharp (like the blue curve above) and occurs faster than I am sampling. There is an example below.

For the saturation, a cheap trick would be to take the left-most and right-most data points. However, I'm not entirely sure how to go about that even though it sounds pretty simple. I could find the minimum and maximum x-coordinates using something like

Min[Table[pairs[[x, 1]], {x, 1, Length[pairs]}]]

but then I don't know how to get the corresponding y-coordinate. Maybe the approach should be very different...

The overall issue is I will be analyzing hundreds of data sets like this and it is completely impractical to try to find these values by hand.

Here is an example of what some of my real data looks like:

Edit: Here is a sample data set in a csv file. I use the following code in Mathematica to import the csv and make the data usable:

imp = Import["sampledata.csv", "Table"];
coilcurrents = Table[.8*imp[[x, 1]], {x, 1, Length[imp]}];
resistances = Table[0.05/imp[[x, 2]], {x, 1, Length[imp]}];
pairs = Table[{1810*coilcurrents[[x]], resistances[[x]]}, {x, 1,
Length[coilcurrents]}];
plot = ListLinePlot[pairs, Frame -> True,
ImageSize -> {575., Automatic}]

The first few lines convert the raw data to meaningful units and can be ignored for all intents and purposes. If anyone is feeling ambitious enough to try something out with this data, be my guest! Added in response to a suggestion.

Final Edit: I've accepted an answer. Here is what my code looks like (the appearance of "25" everywhere is because I have Python spit out my Mathematica code and loop over a bunch of iterations):

imp25 = Import["RHLoopArrayV(25).csv", "Table"];
coilcurrents25 = Table[.8*imp25[[x, 1]], {x, 1, Length[imp25]}];
resistances25 = Table[0.05/imp25[[x, 2]], {x, 1, Length[imp25]}];
pairs25 = Table[{1810*coilcurrents25[[x]], resistances25[[x]]}, {x, 1,
Length[coilcurrents25]}];
pairsordered25 = Sort[pairs25];
pairpairs25 = Partition[pairsordered25, 2, 1];
thresh25 = (Max[pairs25[[All, 2]]] - Min[pairs25[[All, 2]]])/2.
distantcousins25 =
Select[pairpairs25, Abs[#[[1, 2]] - #[[2, 2]]] > thresh25 &];
jumps25 = {distantcousins25[[1, 1, 1]], distantcousins25[[-1, -1, 1]]}
rsat25 = {First[SortBy[pairs25, First]][[2]],
Last[SortBy[pairs25, First]][[2]]}
hmin25 = First[SortBy[pairs12, First]][[1]];
hmax25 = Last[SortBy[pairs12, First]][[1]];
plot25simp =
ListLinePlot[pairs25, Frame -> True,
ImageSize -> {575., Automatic}];
out25=Show[plot25simp,
Graphics[{Thick, Green,
Line[{{jumps25[[1]], rsat25[[1]]}, {jumps25[[1]], rsat25[[2]]}}]}],
Graphics[{Thick, Green,
Line[{{jumps25[[2]], rsat25[[1]]}, {jumps25[[2]], rsat25[[2]]}}]}],
Graphics[{Thick, Green,
Line[{{hmin25, rsat25[[1]]}, {jumps25[[2]], rsat25[[1]]}}]}],
Graphics[{Thick, Green,
Line[{{jumps25[[1]], rsat25[[2]]}, {hmax25, rsat25[[2]]}}]}],FrameLabel -> {"", "", Text[Style["Device 25", Large]], ""}]

This gives me a plot of the real data as well as a representation of the values found by Mathematica.

-
Did you consider fitting a function through your data and then extract the critical values? –  VLC Nov 26 '12 at 18:49
@VLC I considered it but this data does not represent a "function" as it is multi-valued. –  skratch Nov 26 '12 at 18:56
Sort the pairs by x value. Your hysteresis range of interest is where you find neighboring pairs with wildly different y values. To obtain those you could, say partition into consecutive pairs of pairs, then select from those only the ones whose magnitude difference in y values exceeds some given threshold. If you provide sample raw data somebody might put this or another approach into code. –  Daniel Lichtblau Nov 26 '12 at 19:13
@DanielLichtblau I have added raw data as you requested. The data is by nature ordered sequentially by x value starting from zero. The x values end up repeating because the measurement requires that we end up back at the starting x value. I was thinking of something similar to what you're describing but I don't immediately see how to implement it. –  skratch Nov 26 '12 at 19:49
I don't have enough time to flesh this out or improve its efficiency, but I'd consider fitting straight lines to the upper and lower curves shifting the included x-ranges until you maximize a goodness of fit measure. Then, analyze your data relative to those lines. Also, I'd separate your data based on whether your increasing or decreasing the mag field for readability. @whuber may be able to comment on how to do it efficiently. –  rcollyer Nov 26 '12 at 19:59

This is what I had meant. I start with the setup provided in the post and show code from there on.

pairsordered = Sort[pairs];
pairpairs = Partition[pairsordered, 2, 1];
thresh = (Max[pairs[[All, 2]]] - Min[pairs[[All, 2]]])/2.

(* Out[72]= 43.2141 *)

distantcousins =
Select[pairpairs, Abs[#[[1, 2]] - #[[2, 2]]] > thresh &];
jumps = {distantcousins[[1, 1, 1]],
distantcousins[[-1, -1, 1]]}

(* Out[79]= {-114.373, 12.2589} *)

So the hopping points are around -114.4 and 12.3.

-
This looks pretty good... I'll play around with this with a few other data sets and see if there are any hiccups. Thank you! –  skratch Nov 26 '12 at 22:42
For some particularly ugly data this sometimes gets the wrong hopping points but it works overall. Thank you very much for the answer! I may end up playing around a little with the threshold definition to get it to work differently. –  skratch Nov 28 '12 at 16:17

As you haven't provided data let's create some first.

pairs = Table[{RandomReal[], RandomReal[]}, {16}];

Sort them by the x-value and extract the first and the last which will give you the minimum and the maximum for X and Y.

SortBy[pairs, First] // First
(* {0.0997112,0.491322} *)

SortBy[pairs, First] // Last
(* {0.984322,0.653176} *)

Now sample a step function and find the x,y value where it is steepest which should get you close to the point H_c

errorf = Table[{x + 0.2, Erf[x]}, {x, -2, 2, 0.2}];

Differenciate the values in errorf and find the maximum which is in this case the last value pair. You can pick those out by applying Last to the resulting list.

SortBy[Most[
Transpose@{errorf[[All, 1]],
RotateLeft[errorf[[All, 2]]] - errorf[[All, 2]]}], Last]
-
Thanks for the response. Sorting my real data by x value may be problematic in the "hysteretic" region where there are two very different y values. The data comes in naturally ordered as subsequent pairs, so I should be able to skip this step anyways. I have added a link to a csv file of real raw data as an edit at the bottom of the original post. I will take a look at implementing the rest of your suggestion when I get a chance; I'm working on some other things for the moment. –  skratch Nov 26 '12 at 21:11
The numerical sorting to find the "saturated" y values works great so far! The second half of your solution is extremely close to working well, but it doesn't properly deal with the fact that there are two coercive values. It found the negative value with a good deal of accuracy. I may work on splitting my "pairs" variable into one for increasing values and one for decreasing values. –  skratch Nov 26 '12 at 21:29
@skratch look at Daniel's solution that might do what you want in a concise way. –  Matariki Nov 26 '12 at 21:33

You can also use Split:

Split[Sort[pairs], Abs[#1[[2]] - #2[[2]]] > yourthreshold &];

Example:

pairs = RandomReal[1, {16, 2}]
(* {{0.438745, 0.0762158}, {0.129662, 0.777732}, {0.0921677,  0.749584},
{0.415785, 0.418451}, {0.554928, 0.532807}, {0.301706,  0.921911},
{0.729466, 0.750203}, {0.065516, 0.659042}, {0.161117,  0.991978},
{0.144656, 0.805994}, {0.0755413, 0.827263}, {0.447725, 0.0490336},
{0.982973, 0.981803}, {0.553628, 0.521829}, {0.97901, 0.441129},
{0.142663, 0.365815}} *)

sortedpairs = Sort[pairs];
{{xmin, xmax}, {ymin, ymax}} = {Min@#, Max@#} & /@ Transpose[sortedpairs];
threshold = (ymax - ymin)/2;
splitAtJumps = Split[sortedpairs, Abs[#1[[2]] - #2[[2]]] > threshold &];
jumpPoints = First /@ Select[splitAtJumps, Length[#] >= 2 &];
jumPosXcoords = jumpPoints /. {x_, y_} :> {x, 0};
jumPosYcoords = jumpPoints /. {x_, y_} :> {0, y};

In a picture:

Show[ListPlot[splitAtJumps, Joined -> True,
PlotStyle -> Directive[Opacity[.8], Thickness[.015]]],
ListPlot[Join @@ splitAtJumps, Joined -> False,
PlotStyle -> Directive[Red, PointSize[.03]]],
ListPlot[Join @@ splitAtJumps, Joined -> True,
PlotStyle -> Directive[Thick, Dashed, PointSize[.0]]],
ListPlot[jumpPoints , Joined -> False,
PlotStyle -> Directive[Orange, Opacity[.7], PointSize[.05]]],
ListPlot[jumPosXcoords, Joined -> False,
PlotStyle -> Directive[Blue, PointSize[.02]]],
ListPlot[jumPosYcoords, Joined -> False,
PlotStyle -> Directive[Green, PointSize[.02]]],
GridLines -> Transpose[jumpPoints ], ImageSize -> 500,
AxesOrigin -> {0, 0}, PlotRange -> All]

-