0001 function [x] = powerdiff_dke_jd(x1,dx1,x2,dx2,dxmax,y)
0002
0003
0004
0005
0006 if dx2 < dx1
0007 dxc = dx1;
0008 dx1 = dx2;
0009 dx2 = dxc;
0010 inv = 1;
0011 else
0012 inv = 0;
0013 end
0014
0015 deltax = x2 - x1;
0016
0017 N = floor((log(max(dxmax,dx2)) - log(dx1))/log(y)) + 1;
0018
0019 yn = y.^(1:N);
0020 yyn = repmat(yn',[1,N]);
0021 lyn = sum(triu(yyn));
0022
0023 rdxramp = dx2/dx1;
0024 nrdxramp = max([ceil(log(rdxramp)/log(y)) - 1,0]);
0025 if nrdxramp > 1
0026 xramp = (lyn(nrdxramp-1) + 1)*dx2/yn(nrdxramp);
0027 elseif nrdxramp == 1
0028 xramp = dx2/yn(nrdxramp);
0029 else
0030 xramp = 0;
0031 end
0032 solution = 0;
0033 if deltax < max(xramp + dx2/yn(nrdxramp + 1))
0034 if nrdxramp == 1 & deltax >= max(xramp)
0035 if deltax <= dx1*y
0036 solution = 1;
0037 end
0038 end
0039 if nrdxramp > 1 & deltax >= max(xramp)
0040 solution = 1;
0041 end
0042 if solution == 0,
0043 x = [];
0044 return
0045 end
0046 end
0047
0048 dx0 = dxmax;
0049
0050 rdx1 = dx0/dx1;
0051 nrdx1 = floor(log(rdx1)/log(y));
0052
0053 rdx2 = dx0/dx2;
0054 nrdx2 = floor(log(rdx2)/log(y));
0055
0056 if nrdx1 > 0
0057 dx1_max = yn(nrdx1)*dx1;
0058 xramp1 = lyn(nrdx1)*dx1;
0059 else
0060 dx1_max = dx1;
0061 xramp1 = 0;
0062 end
0063 if nrdx2 > 0
0064 dx2_max = yn(nrdx2)*dx2;
0065 xramp2 = lyn(nrdx2)*dx2;
0066 else
0067 dx2_max = dx2;
0068 xramp2 = 0;
0069 end
0070
0071 d = deltax - xramp1 - xramp2;
0072 dx = [];
0073 dx2ramp = [];
0074 dx2_new = dx2/y;
0075 dx1_new = dx1/y;
0076 x = NaN;
0077 while isempty(dx);
0078 if d < 0
0079 if dx2_max >= dx1_max & (nrdx2 > 0)
0080 nrdx2 = nrdx2 - 1;
0081 if nrdx2 > 0
0082 dx2_max = yn(nrdx2)*dx2;
0083 xramp2 = lyn(nrdx2)*dx2;
0084 else
0085 dx2_max = dx2;
0086 xramp2 = 0;
0087 end
0088 elseif nrdx1 > 0
0089 nrdx1 = nrdx1 - 1;
0090 if nrdx1 > 0
0091 dx1_max = yn(nrdx1)*dx1;
0092 xramp1 = lyn(nrdx1)*dx1;
0093 else
0094 dx1_max = dx1;
0095 xramp1 = 0;
0096 end
0097 end
0098 dx0 = min(dx1_max,dx2_max)*y;
0099 else
0100 ndx0 = ceil(d/dx0);
0101 dx0_new = d/(ndx0);
0102 if (dx0_new >= dx1_max) & (dx0_new >= dx2_max | (nrdx2 == 0 & dx0_new >= dx2_new))
0103 dx = [dx1*yn(1:nrdx1),ones(1,ndx0)*dx0_new,dx2ramp,dx2*yn(nrdx2:-1:1)];
0104 break
0105 else
0106 if nrdx2 == 0
0107 if dx1_max < dx2_new
0108 dx2ramp = [dx2_new, dx2ramp];
0109 xramp2 = xramp2 + dx2_new;
0110 dx2_new = dx2_new/y;
0111 elseif nrdx1 > 0
0112 nrdx1 = nrdx1 - 1;
0113 if nrdx1 > 0
0114 dx1_max = yn(nrdx1)*dx1;
0115 xramp1 = lyn(nrdx1)*dx1;
0116 else
0117 dx1_max = dx1;
0118 xramp1 = 0;
0119 end
0120 dx0 = dx1_max*y;
0121 else
0122 if dx0_new > dx1_new
0123 if dx0_new > dx2_new
0124 dx1_max = dx1_new;
0125 elseif dx1_new < dx2_new
0126 dx2ramp = [dx2_new, dx2ramp];
0127 xramp2 = xramp2 + dx2_new;
0128 dx2_new = dx2_new/y;
0129 end
0130 else
0131 warning('empty output from powerdiff');
0132 x = [];
0133 break
0134 end
0135 end
0136 else
0137 if dx1_max >= dx2_max
0138 nrdx1 = nrdx1 - 1;
0139 else
0140 nrdx2 = nrdx2 - 1;
0141 end
0142 if nrdx1 > 0
0143 dx1_max = yn(nrdx1)*dx1;
0144 xramp1 = lyn(nrdx1)*dx1;
0145 else
0146 dx1_max = dx1;
0147 xramp1 = 0;
0148 end
0149 if nrdx2 > 0
0150 dx2_max = yn(nrdx2)*dx2;
0151 xramp2 = lyn(nrdx2)*dx2;
0152 else
0153 dx2_max = dx2;
0154 xramp2 = 0;
0155 end
0156 dx0 = min(dx1_max,dx2_max)*y;
0157 end
0158 end
0159 end
0160 d = deltax - xramp1 - xramp2;
0161 end
0162 if isempty(x)
0163 return
0164 end
0165 if inv == 0;
0166 x = cumsum([x1,dx]);
0167 else
0168 x = cumsum([x1,fliplr(dx)]);
0169 end
0170
0171
0172
0173
0174
0175
0176