# HG changeset patch # User A.M. Thurnherr # Date 1331282854 0 # Node ID 1189eebb260b237e834345c10aea1203a3a0f5cf # Parent 24b3112dc4acdac7c62d7bcf81907fa16270785c . diff -r 24b3112dc4ac -r 1189eebb260b HISTORY --- a/HISTORY Wed Apr 27 11:07:07 2011 -0400 +++ b/HISTORY Fri Mar 09 08:47:34 2012 +0000 @@ -1,9 +1,9 @@ ====================================================================== H I S T O R Y doc: Tue Aug 16 11:41:46 2005 - dlm: Wed Apr 27 11:06:33 2011 + dlm: Fri Jan 6 11:59:11 2012 (c) 2005 A.M. Thurnherr - uE-Info: 164 68 NIL 0 0 72 3 2 4 NIL ofnI + uE-Info: 185 59 NIL 0 0 72 3 2 4 NIL ofnI ====================================================================== NB: CHANGE VERSION IN [default.m] BEFORE UPLOADING @@ -165,3 +165,22 @@ Apr 27, 2011: - updated version in [default.m] to IX_7 - publised on Mercurial server + +IX_8beta: + Jul 2011: + - added new plot options to [geterr.m] + - adapted [getinv.m] to new [geterr.m] + - disabled bin-remapping in [loadrdi.m] + Dec 30, 2011: + - implemented user-provided bug fix in [fixcompass.m] + +IX_8gamma: + Jan 6, 2012: + - removed lesqchol() from getinv.m because it does not + deal well with (near-)singular matrices; while this is + somewhat iffy in the sense that the new code can produce trash + where the old code produced some mix of nans and zeros; + however, if the data are really trash they will cause very + large velocities that are removed from the solutions, i.e. it + is hoped that this change does not confuse things + - added imagnan.m diff -r 24b3112dc4ac -r 1189eebb260b default.m --- a/default.m Wed Apr 27 11:07:07 2011 -0400 +++ b/default.m Fri Mar 09 08:47:34 2012 +0000 @@ -1,9 +1,9 @@ %====================================================================== % D E F A U L T . M % doc: Sat Jun 26 06:10:09 2004 -% dlm: Wed Apr 27 11:05:37 2011 +% dlm: Fri Jan 6 10:58:10 2012 % (c) 2004 ladcp@ -% uE-Info: 28 45 NIL 0 0 72 0 2 4 NIL ofnI +% uE-Info: 28 50 NIL 0 0 72 0 2 4 NIL ofnI %====================================================================== % CHANGES BY ANT: @@ -25,7 +25,7 @@ % the data % structure ps.??? contains parameter for the solution % structure att.??? contains attributes -p.software='LDEO LADCP software: Version IX_7'; +p.software='LDEO LADCP software: Version IX_8gamma'; % file names % f.ladcpdo is the ONLY required input diff -r 24b3112dc4ac -r 1189eebb260b fixcompass.m --- a/fixcompass.m Wed Apr 27 11:07:07 2011 -0400 +++ b/fixcompass.m Fri Mar 09 08:47:34 2012 +0000 @@ -1,3 +1,11 @@ +%====================================================================== +% F I X C O M P A S S . M +% doc: Fri Dec 30 14:59:27 2011 +% dlm: Sat Dec 31 12:32:47 2011 +% (c) 2011 A.M. Thurnherr +% uE-Info: 31 92 NIL 0 0 72 0 2 4 NIL ofnI +%====================================================================== + function [d,p]=fixcompass(d,p) % % fix compass issues @@ -17,6 +25,11 @@ % M. Visbeck LDEO 2004 % +% Changes by ant: +% Dec 30, 2011: - BUG: p.fix_compass = 2 and = 3 did not work; bug report and fix +% provided by Roman Tarakanov today +% bug fix verified with PALE data (old code produces garbage results) + disp('FIXCOMPASS adjust compass') % use tilt sensors to figure out allingment between instruments @@ -43,7 +56,7 @@ if p.fix_compass==3 disp(' use up looking compass on down looking ADCP') - rotv(1,:)=d.hdg(1,:)-d.hdg(2,:)+rot(1); + rotv(1,:)=d.hdg(2,:)-d.hdg(1,:)+rot(1); p.rotup2down=0; end @@ -53,7 +66,7 @@ p.hdg_offset_used(2)=p.hdg_offset_used(2)+rot(2); if p.fix_compass==2 disp(' use down looking compass on up looking ADCP') - rotv(2,:)=d.hdg(2,:)-d.hdg(1,:)+rot(2); + rotv(2,:)=d.hdg(1,:)-d.hdg(2,:)+rot(2); p.rotup2down=0; end end diff -r 24b3112dc4ac -r 1189eebb260b geterr.m --- a/geterr.m Wed Apr 27 11:07:07 2011 -0400 +++ b/geterr.m Fri Mar 09 08:47:34 2012 +0000 @@ -1,9 +1,9 @@ %====================================================================== % G E T E R R . M % doc: Wed Jun 30 23:24:51 2004 -% dlm: Fri Mar 5 15:48:05 2010 +% dlm: Wed Jul 6 20:26:14 2011 % (c) 2004 ladcp@ -% uE-Info: 91 17 NIL 0 0 72 2 2 8 NIL ofnI +% uE-Info: 19 35 NIL 0 0 72 2 2 8 NIL ofnI %====================================================================== % MODIFICATIONS BY ANT: @@ -13,14 +13,22 @@ % Jul 16, 2004: - added global variable skip_figure_3 to workaround % linux matlab bug % Oct 7, 2008: - extensively modified procfig 3 for version IX_6 +% Jun 29, 2011: - removed skp_figure_3 +% - added ps.fig3_colormap, ps.fig3_err_y_axis, ps.fig3_avgerr +% Jun 30, 2011: - fixed fig.3 middle column plot title for median plot +% Jul 6, 2001: - fixed plot title -function l=geterr(dr,d,iplot) +function l=geterr(ps,dr,d,iplot) % function l=geterr(dr,d,iplot) % returns predicitons of U_ocean and % U_ctd on the raw data grid % % CTD velocity -if nargin<3, iplot=1; end +if nargin<4, iplot=1; end + +ps=setdefv(ps,'fig3_colormap',2); % 1: jet 2: polar +ps=setdefv(ps,'fig3_err_y_axis',2); % 1: bin# 2: depth +ps=setdefv(ps,'fig3_avgerr',2'); % 1: mean 2: median tim=dr.tim; tim(1)=-1e30; @@ -102,7 +110,9 @@ iiz=ceil(iz(ii(j))); iit=itm(ii(j)); l.u_oce(iiz,iit)=d.ru(ii(j))-l.ru_ctd(ii(j)); + l.u_err(iiz,iit)=d.ru(ii(j))-l.ru_oce(ii(j))-l.ru_ctd(ii(j)); l.v_oce(iiz,iit)=d.rv(ii(j))-l.rv_ctd(ii(j)); + l.v_err(iiz,iit)=d.rv(ii(j))-l.rv_oce(ii(j))-l.rv_ctd(ii(j)); l.w_oce(iiz,iit)=d.rw(ii(j))-l.rw_ctd(ii(j)); l.w_oce_z(iiz,iit)=d.rw(ii(j))-l.rw_ctd_z(ii(j)); if existf(l,'rw_ctd_p') @@ -120,7 +130,9 @@ end else % uplooker and downlooker bin sizes are equal l.u_oce=full(sparse(ceil(iz(ii)),itm(ii),d.ru(ii)-l.ru_ctd(ii))); + l.u_err=full(sparse(ceil(iz(ii)),itm(ii),d.ru(ii)-l.ru_oce(ii)-l.ru_ctd(ii))); l.v_oce=full(sparse(ceil(iz(ii)),itm(ii),d.rv(ii)-l.rv_ctd(ii))); + l.v_err=full(sparse(ceil(iz(ii)),itm(ii),d.rv(ii)-l.rv_oce(ii)-l.rv_ctd(ii))); l.w_oce=full(sparse(ceil(iz(ii)),itm(ii),d.rw(ii)-l.rw_ctd(ii))); l.w_oce_z=full(sparse(ceil(iz(ii)),itm(ii),d.rw(ii)-l.rw_ctd_z(ii))); if existf(l,'rw_ctd_p') @@ -177,41 +189,69 @@ d.ru(ii)=nan; d.rv(ii)=nan; -global skip_figure_3; -if isempty(skip_figure_3) - figure(3) clf orient landscape - colormap(polarmap(21)); - - subplot(231) + + if ps.fig3_colormap == 2 + colormap(polarmap(21)); + else + col=jet(128); + col=([[1 1 1]; col]); + colormap(col) + end + ib=1:size(l.ru_err,1); ib=ib-length(d.izu); - tmp = l.ru_err; tmp(isnan(tmp)) = 0; - pcolorn(l.itv2,-ib,tmp), shading flat + + subplot(231) + if ps.fig3_err_y_axis == 2 + if ps.fig3_colormap == 2 + tmp = l.u_err; tmp(isnan(tmp)) = 0; + pcolorn(l.itv,-l.z_oce,tmp), shading flat + else + pcolorn(l.itv,-l.z_oce,l.u_err), shading flat + end + ylabel('Depth [m]'); + else + if ps.fig3_colormap == 2 + tmp = l.ru_err; tmp(isnan(tmp)) = 0; + pcolorn(l.itv2,-ib,tmp), shading flat + else + pcolorn(l.itv2,-ib,l.ru_err), shading flat + end + ylabel('Bin #'); + end fac=meannan(l.u_oce_s); fac=max([fac,1e-2]); caxis([-3 3]*fac) colorbar xlabel('Super Ensemble #'); - ylabel('Bin #'); title(sprintf('U-err std: %.03f',meannan(stdnan(l.ru_err')))) subplot(232) - plot(meannan(l.ru_err')',-ib) + if ps.fig3_avgerr == 2 + plot(medianan(l.ru_err')',-ib) + title('median(U-err)') + else + plot(meannan(l.ru_err')',-ib) + title('mean(U-err)') + end set(gca,'XLim',[-0.05 0.05]); set(gca,'Ylim',[-ib(end) -ib(1)]); set(gca,'Xtick',[-0.04:0.02:0.04]); grid xlabel('Residual [m/s]'); ylabel('Bin #'); - title('mean(U-err)') subplot(233) - tmp = l.u_oce; tmp(isnan(tmp)) = 0; - pcolorn(l.itv,-l.z_oce,tmp), shading flat + if ps.fig3_colormap == 2 + tmp = l.u_oce; tmp(isnan(tmp)) = 0; + pcolorn(l.itv,-l.z_oce,tmp), shading flat + else + pcolorn(l.itv,-l.z_oce,l.u_oce), shading flat + end ca = caxis; if abs(ca(1)) > abs(ca(2)) caxis([-abs(ca(1)) abs(ca(1))]); @@ -231,8 +271,23 @@ title('U_{oce}') subplot(234) - tmp = l.rv_err; tmp(isnan(tmp)) = 0; - pcolorn(l.itv2,-ib,tmp), shading flat + if ps.fig3_err_y_axis == 2 + if ps.fig3_colormap == 2 + tmp = l.v_err; tmp(isnan(tmp)) = 0; + pcolorn(l.itv,-l.z_oce,tmp), shading flat + else + pcolorn(l.itv,-l.z_oce,l.v_err), shading flat + end + ylabel('Depth [m]'); + else + if ps.fig3_colormap == 2 + tmp = l.rv_err; tmp(isnan(tmp)) = 0; + pcolorn(l.itv2,-ib,tmp), shading flat + else + pcolorn(l.itv2,-ib,l.rv_err), shading flat + end + ylabel('Bin #'); + end fac=meannan(l.v_oce_s); fac=max([fac,1e-2]); caxis([-3 3]*fac) @@ -242,18 +297,27 @@ title(sprintf('V-err std: %.03f',meannan(stdnan(l.rv_err')))) subplot(235) - plot(meannan(l.rv_err')',-ib) + if ps.fig3_avgerr == 2 + plot(medianan(l.rv_err')',-ib) + title('median(V-err)') + else + plot(meannan(l.rv_err')',-ib) + title('mean(V-err)') + end set(gca,'XLim',[-0.05 0.05]); set(gca,'Ylim',[-ib(end) -ib(1)]); set(gca,'Xtick',[-0.04:0.02:0.04]); grid xlabel('Residual [m/s]'); ylabel('Bin #'); - title('mean(V-err)') subplot(236) - tmp = l.v_oce; tmp(isnan(tmp)) = 0; - pcolorn(l.itv,-l.z_oce,tmp), shading flat + if ps.fig3_colormap == 2 + tmp = l.v_oce; tmp(isnan(tmp)) = 0; + pcolorn(l.itv,-l.z_oce,tmp), shading flat + else + pcolorn(l.itv,-l.z_oce,l.v_oce), shading flat + end ca = caxis; if abs(ca(1)) > abs(ca(2)) caxis([-abs(ca(1)) abs(ca(1))]); @@ -270,10 +334,9 @@ colorbar xlabel('Ensemble #'); ylabel('Depth [m]'); - title('U_{oce}') + title('V_{oce}') streamer([dr.name,' Figure 3']); -end % of ~skip_figure_3 % reset colormap figure(11) diff -r 24b3112dc4ac -r 1189eebb260b getinv.m --- a/getinv.m Wed Apr 27 11:07:07 2011 -0400 +++ b/getinv.m Fri Mar 09 08:47:34 2012 +0000 @@ -8,9 +8,9 @@ %====================================================================== % G E T I N V . M % doc: Thu Jun 17 15:36:21 2004 -% dlm: Wed Aug 25 02:04:54 2010 +% dlm: Fri Jan 6 11:57:45 2012 % (c) 2004 ladcp@ -% uE-Info: 988 0 NIL 0 0 72 2 2 8 NIL ofnI +% uE-Info: 31 50 NIL 0 0 72 2 2 8 NIL ofnI %====================================================================== % CHANGE HISTORY: @@ -25,6 +25,10 @@ % (and corresponding velocities removed) even % on sadcpfac=0, i.e. when SADCP data were only % used for plotting +% Jun 29, 2011: - pass ps to geterr to allow controlling fig. 3 +% Jan 6, 2012: - lesqchol removed because it does not deal with +% nearly singular matrices gracefully +% - BUG: replaced all imag() by new imagnan() if nargin<5, iplot=0; end @@ -364,7 +368,7 @@ if length(ubot)>0 dr.zbot=z(iubot); dr.ubot=real(ubot(:,1)); - dr.vbot=imag(ubot(:,1)); + dr.vbot=imagnan(ubot(:,1)); dr.uerrbot=ubot(:,2); % make length of array unique while length(dr.zbot)==nt | length(dr.zbot)==nz @@ -508,7 +512,7 @@ % save results in output array dr.z=z; dr.u=real(uocean(:,1)); -dr.v=imag(uocean(:,1)); +dr.v=imagnan(uocean(:,1)); if size(uocean,2)>1 dr.uerr=(uocean(:,2)); end @@ -523,12 +527,12 @@ dr.xship=(slon-slon(1))*60*1852*cos(meannan(slat)*pi/180); dr.yship=(slat-slat(1))*60*1852; dr.uship=real(shipvel); - dr.vship=imag(shipvel); + dr.vship=imagnan(shipvel); end dr.zctd=zctd; dr.wctd=-wctd; dr.uctd=-real(uctd(:,1))'; -dr.vctd=-imag(uctd(:,1))'; +dr.vctd=-imagnan(uctd(:,1))'; if size(uctd,2)>1 dr.uctderr=(uctd(:,2))'; end @@ -537,14 +541,14 @@ dt=mean([0,dt;dt,0]); ctdpos=-cumsum(uctd(:,1).*dt').'; dr.xctd=real(ctdpos); -dr.yctd=imag(ctdpos); +dr.yctd=imagnan(ctdpos); figure(7) clf orient tall tim=dr.tim-fix(dr.tim(1)); uctd_drag=real(ctdvel); - vctd_drag=imag(ctdvel); + vctd_drag=imagnan(ctdvel); % plot some of the drag fac results subplot(421) plot(tim,dr.uctd,'-b','linewidth',1.8) @@ -641,7 +645,7 @@ pause(0.01) % compute velcoity error - der=geterr(dr,di,iplot); + der=geterr(ps,dr,di,iplot); if size(uocean,2)>1 % second pass: uerr is set % note that dr.uerr is entirely re-scaled %%##%% dr.uerr=dr.uerr/medianan(dr.uerr)*medianan(sqrt(der.u_oce_s.^2+der.v_oce_s.^2)); @@ -730,7 +734,7 @@ uocean_do(ii)=NaN; end dr.u_do=real(uocean_do(:,1)); - dr.v_do=imag(uocean_do(:,1)); + dr.v_do=imagnan(uocean_do(:,1)); if length(iupc)>5, @@ -764,7 +768,7 @@ uocean_up=uocean_do*NaN; end dr.u_up=real(uocean_up(:,1)); - dr.v_up=imag(uocean_up(:,1)); + dr.v_up=imagnan(uocean_up(:,1)); u_bias = meannan(dr.u_do-dr.u_up); v_bias = meannan(dr.v_do-dr.v_up); @@ -864,8 +868,6 @@ Ab(:,id)=[]; if length(ibot)>1 % solve - %[ubot,uebot]=lesqfit(dbot,Ab); - % ubot=lesqchol(dbot,Ab); [m,me]=lesqfit(dbot,Ab); ubot=[full(m),abs(full(me))]; velerr=min(abs(full(me))); @@ -978,14 +980,14 @@ ds.z_sadcp=zsadcp(iok); ds.u_sadcp=real(dsadcp(iok)); - ds.v_sadcp=imag(dsadcp(iok)); + ds.v_sadcp=imagnan(dsadcp(iok)); ds.uerr_sadcp=verr(iok); end % nargout <= 3 else % sadcpfac <= 0 if nargout>3 ds.z_sadcp=zsadcp; ds.u_sadcp=real(dsadcp); - ds.v_sadcp=imag(dsadcp); + ds.v_sadcp=imagnan(dsadcp); ds.uerr_sadcp=verr; end end @@ -1170,14 +1172,14 @@ function [uocean,uctd,uoceanerr,uctderr]=lainsolv(Aocean,Actd,dladcp,nsolve); %function [uocean,uctd,uoceanerr,uctderr]=lainsolv(Aocean,Actd,dladcp,nsolve); % -% solve LADCP current profiling using an invers technique. +% solve LADCP current profiling using an inverse technique. % dladcp = [Aocean,Actd] * [uocean,uctd] % % nsolve = 0 Cholseky transform % = 1 Moore Penrose Inverse % % Martin Visbeck LDEO 15/12/1998 -% need LESQFIT and LESQCHOL routines +% need LESQFIT routine % if nargin<4, nsolve=0; end @@ -1196,8 +1198,8 @@ [m,me]=lesqfit(d,A); me=full(me); else - disp('Cholesky transform') - m=lesqchol(d,A); + disp('Moore-Penrose inverse w/o errors') + m=lesqfit(d,A); end % split results up @@ -1241,31 +1243,6 @@ Ac=Ac.*sparse(ii,ic,w); %------------------------------------------------------------------- -function [m,dm,c]=lesqchol(d,g) -% function [m,dm,c]=lesqcholw(d,g) - -% fit least squares method to linear problem -% Use Cholesky transform -% -%input parameters: -% d:= data vector ; g:= model matrix -% output parameters: -% m=model factors; dm= model data, c=correlation - -n=length(d); -[i,j]=size(g); -if i~=n; disp(' wrong arguments'),return,end -[r,b] = chol( g.' * g); -if b~=0, m=g(1,:)'+NaN; dm=d+NaN; c=NaN; return, end -y = forwardsub(r.' , g.' * d); -m = backsub(r,y); -if nargout<2, return, end -dm = g * m; -if nargout<3, return, end -co = cov([d,dm]); -c = co(1,2) / sqrt( co(1,1)*co(2,2) ); - -%------------------------------------------------------------------- function [m,me,c,dm,gi]=lesqfit(d,g) % function [m,me,c,dm,gi]=lesqfit(d,g) diff -r 24b3112dc4ac -r 1189eebb260b imagnan.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/imagnan.m Fri Mar 09 08:47:34 2012 +0000 @@ -0,0 +1,14 @@ +%====================================================================== +% I M A G N A N . M +% doc: Fri Jan 6 11:31:23 2012 +% dlm: Fri Jan 6 11:33:24 2012 +% (c) 2012 A.M. Thurnherr +% uE-Info: 10 48 NIL 0 0 72 0 2 8 NIL ofnI +%====================================================================== + +% HISTORY: +% Jan 6, 2012: - created for version IX_8gamma + +function I = imagnan(C) + I = imag(C); + I(find(isnan(real(C)))) = nan; diff -r 24b3112dc4ac -r 1189eebb260b loadrdi.m --- a/loadrdi.m Wed Apr 27 11:07:07 2011 -0400 +++ b/loadrdi.m Fri Mar 09 08:47:34 2012 +0000 @@ -13,9 +13,9 @@ %====================================================================== % L O A D R D I . M % doc: Fri Jun 18 18:21:56 2004 -% dlm: Wed Oct 28 11:23:27 2009 +% dlm: Thu Aug 18 17:20:09 2011 % (c) 2004 ladcp@ -% uE-Info: 1574 0 NIL 0 0 72 2 2 8 NIL ofnI +% uE-Info: 46 76 NIL 0 0 72 2 2 8 NIL ofnI %====================================================================== % CHANGES BY ANT @@ -40,6 +40,10 @@ % Nov 18, 2007: - added code for p.allow_3beam_solutions, p.ignore_beam % Jan 7, 2009: - tightened use of exist() % Oct 28, 2009: - modified p.ignore_beam for dual-head systems +% Jun 27, 2011: - l.blen removed because bin lenght can be different for UL & DL +% - apparently unused z-variable commented out +% Jun 30, 2011: - buggy bin-remapping disabled +% Aug 18, 2011: - added comment to coord-transformation code (gimbal pitch) % p=setdefv(p,'pg_save',[1 2 3 4]); % Default =3 for loadctd_whoi. @@ -587,7 +591,7 @@ % L = UPDOWN('filedown','fileup') reads ADCP raw data from the specified % files. L is a structure array with the following fields: % -% blen: bin length +% blen: bin length %%% ANT: REMOVED 2011/06/28 BECAUSE IT CAN BE DIFFERENT FOR UL/DL % nbin: number of bins % blnk: blank after transmit % dist: distance of bin 1 from transducer @@ -848,7 +852,7 @@ end % distance vectors -z = fd(f.dist) + fd(f.blen)*([1:fd(f.nbin)] - 1); +%%% z = fd(f.dist) + fd(f.blen)*([1:fd(f.nbin)] - 1); % unused?! ANT 2011/06/28 if bmax(1)>0, fd(f.nbin)=bmax(1); end idb=1:fd(f.nbin); @@ -945,7 +949,6 @@ l.npng_d = fd(f.npng); l.nens_u = size(vu,1); l.nens_d = size(vd,1); - l.blen = fd(f.blen); l.nbin = fd(f.nbin); l.blnk = fd(f.blnk); l.dist = fd(f.dist); @@ -1011,11 +1014,10 @@ l.problem = [fliplr(uproblem(iu,:)) dproblem(id,:)]'; l.problemb = problemb(id,:)'; -else +else % single instrument l.zd=[0:(fd(f.nbin)-1)]*fd(f.blen)+fd(f.dist); eval(['l.serial_cpu_d=fd(',f.serial,');']); - l.blen = fd(f.blen); l.npng_d = fd(f.npng); l.nens_d = length(vd(v.tim)); l.nbin = fd(f.nbin); @@ -1566,7 +1568,7 @@ a=setdefv(a,'use_tilt',1); a=setdefv(a,'use_heading',1); -a=setdefv(a,'use_binremap',1); +a=setdefv(a,'use_binremap',0); %%% CODE IS BUGGY! DO NOT USE a=setdefv(a,'beams_up',a.Up); a=setdefv(a,'beamangle',a.Beam_angle); @@ -1616,6 +1618,14 @@ RR=roll.*d2r; KA=sqrt(1.0 - (sin(pitch.*d2r).*sin(roll.*d2r)).^2); PP=asin(sin(pitch.*d2r).*cos(roll.*d2r)./KA); +%% NB: The preceding two lines could be replaced with +%% PP=atan(tan(pitch.*d2r) * cos(roll.*d2r)); +%% which is the expression given by RDI in the coord- +%% trans manual. I have tried this with a single +%% file from DIMES UK2 and the max velocity differences +%% are 1e-13 m/s, i.e. they look consistent with +%% roundoff errors. + HH=head.*d2r; % Step 2 - calculate trig functions and scaling factors @@ -1662,7 +1672,8 @@ % Step 3: correct depth cell index for pitch and roll for i=1:4, if a.use_binremap - J(i)=fix(IB.*SC(i)+0.5); +% J(i)=fix(IB.*SC(i)+0.5); + J(i)=IB; %%% else J(i)=IB; end