Beta version. Published: 20170404. Updated: 20170410.
Update 20191228: corrected formulas in coilcraft_model.ckt for compatibility (avoiding the ** operator).
Update 20210404: Parameters changed to engineering units. Added mouse drag adjustment (try click and drag on the "tweaked" values!).
Coilcraft has a wonderful selection of inductors, transformers, and related stuff. Most of their catalog is supported with rich data, including SPICE models that are based on real measurements. (I'm told they basically set up an intern with a VNA and fixture, and pounded out their whole catalog in a summer. Well, that's probably not true, but it's amusing to imagine.) There's one big downside, though. Their models use Warburg resistors: a pure resistance, proportional to √f.
This is the equivalent circuit used. The dependent elements are RVAR1 and RVAR2. (LVAR is also a dependent element, but the dependency is usually small in the given frequency range, and can be ignored.)
SPICE Transient Analysis does very poorly with a frequencydependent element. LTSpice is able to analyze the model (Coilcraft provides instructions on how to set this up), but the numerical accuracy is poor, and output is very slow to generate. PSpice also supports this model, probably with the same drawbacks. Some ignore the elements entirely: Multisim 11.0 works in AC Analysis, but in Transient Analysis, it's an open circuit. Altium doesn't appear to support arbitrary transfer functions at all.
Quick Start Instructions:
1. Copy the model parameters from the datasheet, and paste them in the text box below.
2. Check the unit multipliers in the datasheet. Compare them to the multipliers in the format box below. If any do not match, correct the numbers to the required multiplier. (Do not write unit multipliers in the box.)
Two text formats are understood:
<Part Number> <R1 (Ω)> <R2 (Ω)> <C (pF)> <L (nH)> <k> <Fmax (MHz)>
Seven parameters: this is an RF Inductor model. Note the order of parameters, no F_{min}, and L is in nH.
<Part Number> <Fmin (MHz)> <Fmax (MHz)> <R1 (Ω)> <R2 (Ω)> <C (pF)> <k1> <k2> <k3  L (μH)> [k4 [k5]]
Nine or more parameters: this is a Power Inductor model. Note the order of parameters, k3 (or L) is in μH, and k4 and k5 are optional. (Additional parameters are ignored.)
7. This code calls out coilcraft_model.ckt. Make sure it's available on the simulation path, or is copied into the same file.
As far as I know, Coilcraft provides five varieties of SPICE models. Three use the Warburg resistor equivalent model, and are supported by this tool. Examples are listed below, showing one part from each supported format. The remaing two types, which are currently unsupported in this tool, are the transmission line equivalent model (used for RF spring inductors), and the oddball SLC7530D, which is a dual version (the circuit is the same Warburg equivalent model, doubled up; they can still be used, with some modification, assuming you don't depend on coupling between the two windings). There's also Sparameter data, but that's much more difficult to work with.
RF Air Core  

Part Number  R1  R2  C  L, nH  k  Fmax, MHz  
1206CS221  25  0.5  0.082  223  3.15E04  1300  
Power, and Cored RF  
Part Number  F_{min}, MHz  F_{max}, MHz  R_{1}, Ω  R_{2}, Ω  C, pF  k1  k2  k3 or L, uH  k4  k5 
DO3308P104  4  8  0.076  28.6  2.65  2.15E07  16.9  94.0  
XAL4020222  0.1  100  100  0.0352  4.7  1.6E04  0.48  2.2  0.01  8.0E06 
Follow the numbered steps above. If an input is entered wrong, a warning will pop up when the button is pressed. Numeric inputs cannot have letters in them (besides scientific notation). Check the inputs for the correct unit multiplier, and shift the numbers as needed. Some parameters will fail softly: an error in Fmin, k1, k2, k4 or k5 will result in using the default value instead. The defaults are: Fmin = Fmax × 10^{4}, k1 = 1e9, k2 = 1e9, k4 = 0, and k5 = 1e9. This happens most often when pasting in a model that has missing parameters, and these defaults should be fine.
When pressing the Update button, the Parsed Values column will only be updated when no errors are found. The Tweaked Values column is only updated when the Copy button is pressed, so you can work with two different models at once. (Keep this in mind if you are working to match models.)
The Copy button automatically refreshes the plot and model parameter output, as does pressing any selection button, or the Refresh button. The output is only updated if no errors are found in the Tweaked values.
The vertical axis is autoranging, clamped to some basic bounds (1e±37) to keep dividebynearzero results from blowing up too badly. In addition, Q is clamped to a minimum of 1, because below 1, the impedance is more resistor than inductor. The axes are always loglog (this may be made selectable later).
A Warburg element arises in the equivalent circuit of a diffusion process. Examples include ionic diffusion (a battery's terminal impedance), semiconductor charge diffusion (which rules the behavior of diodes and BJTs), or the diffusion of currents into conductors (here, skin effect and core loss).
Thing is, they use a frequencydependent resistor. This has a nonphysical phase shift: no real twoterminal device can exhibit impedance variation versus frequency, without a corresponding phase shift. (The proper definition of a Warburg element has the phase shift at 45°: half reactive, half resistive.) Also, an ideal Warburg element has an indefinite impulse response. This works great in frequencydomain calculations (AC Analysis), but makes most simulators hopelessly fall over in the time domain (Transient Analysis).
The problem is: to perform a Transient Analysis of a frequencydependent resistor, the simulator has to calculate the impulse response of the network (its inverse Laplace transform). To calculate a timestep, the response is convolved with the running history of the transient simulation. Well, because the impulse response is ideally infinite, the simulator has to make some compromise, limiting the length of that response. Thus, for short enough time scales (high frequencies), it will be infinite (or more likely, RSHUNT, or 1 / GMIN), and at very long time scales (very low frequencies), it will eventually fall to zero. Second, it's accidentally quadratic (or perhaps not accidentally, since that's simply the price to pay for convolution). The convolution must be done over the full history of timesteps, for each timestep. So even if you get a simulation that works (doesn't crash out with "timestep too small", yields reasonable looking output), it will necessarily go very, very slowly.
So what can we substitute with? When a diffusion element is needed, in circuit, it is either approximated with a finite network, or computed with an ADCDSPDAC system. An example of this appears in The Art of Electronics, 3rd ed., p.559, where an RC network is used to approximate a 1/√f response. The design is simple: a parallel set of R+C, or a series chain of R∥C, where the R and C values are distributed geometrically. The poles are real, spaced evenly, and the impedance has equal parts resistance and reactance for all frequencies in the range. For the inverse function (Z rising with f), use L instead of C.
A more detailed example can be found in: A comprehensive study on different approximation methods of fractional order system, Iqbal and Shekh, International Research Journal of Engineering and Technology (IRJET) Vol 03 Iss 08 (Aug 2016), including discussion of analysis, and approximation methods for other exponents.
Since we need an impedance rising with frequency, we'll choose this network. This approximates an actual Warburg element (with the natural phase shift), over a limited frequency range. The transfer function consists of N poles and zeroes; we can space the polezero pairs so that the impedance stays very close to the ideal √f curve, within the desired range. Outside of that range, the impedance becomes asymptotic: inductive at low frequencies, resistive at high frequencies. One more (unpaired) resistor, or inductor, can be added in series to change either asymptote, if an inductivetoinductive or resistivetoresistive response is desired. We'll take advantage of this, in a little bit.
An approximated Warburg element takes care of RVAR1 and RVAR2. Unaccounted for is LVAR. The original model has a logarithmic correction versus frequency, usually slightly larger at very low frequencies (but trending towards infinity at DC, which is quite peculiar), and slightly smaller at high frequencies. On the other hand, some use an opposite coefficient (more inductance at higher frequencies), which is also very peculiar. Either way, resistance is not factored into this adjustment (so violates Kramers–Kronig relations). For most models (L decreasing with frequency), the decrease is moreorless compensated with the Warburg elements. For the others, L and Q will be underestimated, and the curves will match poorly.
The most common disparity is in the SRF (selfresonant frequency), which is usually higher than the original model. This is because the reactance of the Warburg element acts in parallel with the inductor, reducing total inductance. Likewise, the apparent inductance rises at low frequency, due to RVAR1: but this is not a mistake, simply a quirk of calculating the inductance when the impedance is largely resistive (Q < 1). (To be exact, it's because RVAR1's lumped equivalent is made of inductors and resistors, so that those inductors act in series with the main inductor at very low frequencies.)
A few words about the exact network we'll be substituting with. While a geometric distribution of R and L does a good job, it doesn't do the best job, particularly if an alternative asymptote is chosen. I used the engineer's old friend, Excel, to optimize the values for a 4decade, 5inductor network, normalized for symmetry around 1Hz and 1Ω. This took some preparation and adjustment. In the end, these values were found:


R_{0} is the lone resistor, added in series, to make the LF asymptote resistive.
The impedance is plotted below:
The ideal impedance reference is shown in blue. And the error, as a percentage of impedance:
As you can see, an unusually good fit, for the basic solver that comes with Excel. It takes much refinement to be able to walk a complex, multidimensional, polynomial system towards such a symmetrical result!
Encapsulated in SPICE, and normalized for the LF asymptote, we have rdiff.ckt:
* Diffusion Element * By Tim Williams, 20170330 * * Features: *  Four decade frequency range. *  Resistive LF and HF asymptotes. *  Normalized to LF asymptote (1 ohm, 1 Hz). *  HF asymptote: 175.9534 ohm, 30.95962k Hz. *  Center: 13.26588 ohm, 175.9836 Hz (default parameters * normalize the middle of the range to 1 ohm, 1Hz). *  Accuracy: +/ 0.93% RMS, +/ 1.3% peak, * from (1.33 ohm, 1.76 Hz) to (133 ohms, 17.6k Hz). * .SUBCKT RDIFF 1 7 PARAMS: Resistance=75.38135m Frequency=5.683323m R0 1 2 { 1.0 * Resistance } R1 2 3 { 1.436193 * Resistance } R2 3 4 { 3.936277 * Resistance } R3 4 5 { 10.61526 * Resistance } R4 5 6 { 29.09347 * Resistance } R5 6 7 { 129.8723 * Resistance } L1 2 3 { 41.00917m * Resistance / Frequency } L2 3 4 { 15.51648m * Resistance / Frequency } L3 4 5 { 5.876198m * Resistance / Frequency } L4 5 6 { 2.259226m * Resistance / Frequency } L5 6 7 { 1.2241m * Resistance / Frequency } .ENDS