Index: /issm/trunk-jpl/src/m/contrib/oasis/carter.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/oasis/carter.m (revision 13011) +++ /issm/trunk-jpl/src/m/contrib/oasis/carter.m (revision 13011) @@ -0,0 +1,28 @@ +function attenuation=carter(temperature) +%CARTER - attenuation as a function of temperature +% +% TWO WAY - Attenuation (in dB/m) as a function of temperature (K) +% From Carter at al. 2007 (Radar-based subglacial lake classification in Antarctica) +% Figure 4 +% +% Usage: +% attenuation=carter(temperature) + +if(temperature<0) + error('input temperature should be in Kelvin (positive)'); +end +T=temperature-273.15; + +Temp=[-50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0]'; +A=[0.0015 0.002 0.003 0.0042 0.0055 0.0083 0.012 0.0175 0.026 0.038 0.055]'; + +%Now, do a cubic fit between Temp and B: +[cfun,gof,output]=fit(Temp,A,'cubicspline'); +%breaks=cfun.p.breaks; +%coeff=cfun.p.coefs; + +%Calculate attenuation +attenuation=cfun(T); + +%Make it a 2 way attenuation +attenuation=2*attenuation; Index: /issm/trunk-jpl/src/m/materials/paterson.m =================================================================== --- /issm/trunk-jpl/src/m/materials/paterson.m (revision 13011) +++ /issm/trunk-jpl/src/m/materials/paterson.m (revision 13011) @@ -0,0 +1,42 @@ +function rigidity=paterson(temperature) +%PATERSON - figure out the rigidity of ice for a given temperature +% +% rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). +% temperature is in Kelvin degrees +% +% Usage: +% rigidity=paterson(temperature) + +if(temperature<0) + error('input temperature should be in Kelvin (positive)'); +end +T=temperature-273.15; + +%The routine below is equivalent to: + +% n=3; T=temperature-273; +% %From paterson, +% Temp=[0;-2;-5;-10;-15;-20;-25;-30;-35;-40;-45;-50]; +% A=[6.8*10^-15;2.4*10^-15;1.6*10^-15;4.9*10^-16;2.9*10^-16;1.7*10^-16;9.4* +% 10^-17;5.1*10^-17;2.7*10^-17;1.4*10^-17;7.3*10^-18;3.6*10^-18];;%s-1(kPa-3) +% %Convert into rigidity B +% B=A.^(-1/n)*10^3; %s^(1/3)Pa +% %Now, do a cubic fit between Temp and B: +% fittedmodel=fit(Temp,B,'cubicspline'); +% rigidity=fittedmodel(temperature); + +rigidity=zeros(length(T),1); +pos1=find(T<=-45); rigidity(pos1)=10^8*(-0.000292866376675*(T(pos1)+50).^3+ 0.011672640664130*(T(pos1)+50).^2 -0.325004442485481*(T(pos1)+50)+ 6.524779401948101); +pos2=find(-45<=T & T<-40); rigidity(pos2)=10^8*(-0.000292866376675*(T(pos2)+45).^3+ 0.007279645014004*(T(pos2)+45).^2 -0.230243014094813*(T(pos2)+45)+ 5.154964909039554); +pos3=find(-40<=T & T<-35); rigidity(pos3)=10^8*(0.000072737147457*(T(pos3)+40).^3+ 0.002886649363879*(T(pos3)+40).^2 -0.179411542205399*(T(pos3)+40)+ 4.149132666831214); +pos4=find(-35<=T & T<-30); rigidity(pos4)=10^8*(-0.000086144770023*(T(pos4)+35).^3+ 0.003977706575736*(T(pos4)+35).^2 -0.145089762507325*(T(pos4)+35)+ 3.333333333333331); +pos5=find(-30<=T & T<-25); rigidity(pos5)=10^8*(-0.000043984685769*(T(pos5)+30).^3+ 0.002685535025386*(T(pos5)+30).^2 -0.111773554501713*(T(pos5)+30)+ 2.696559088937191); +pos6=find(-25<=T & T<-20); rigidity(pos6)=10^8*(-0.000029799523463*(T(pos6)+25).^3+ 0.002025764738854*(T(pos6)+25).^2 -0.088217055680511*(T(pos6)+25)+ 2.199331606342181); +pos7=find(-20<=T & T<-15); rigidity(pos7)=10^8*(0.000136920904777*(T(pos7)+20).^3+ 0.001578771886910*(T(pos7)+20).^2 -0.070194372551690*(T(pos7)+20)+ 1.805165505978111); +pos8=find(-15<=T & T<-10); rigidity(pos8)=10^8*(-0.000899763781026*(T(pos8)+15).^3+ 0.003632585458564*(T(pos8)+15).^2 -0.044137585824322*(T(pos8)+15)+ 1.510778053489523); +pos9=find(-10<=T & T<-5); rigidity(pos9)=10^8*(0.001676964325070*(T(pos9)+10).^3- 0.009863871256831*(T(pos9)+10).^2 -0.075294014815659*(T(pos9)+10)+ 1.268434288203714); +pos10=find(-5<=T & T<-2); rigidity(pos10)=10^8*(-0.003748937622487*(T(pos10)+5).^3+0.015290593619213*(T(pos10)+5).^2 -0.048160403003748*(T(pos10)+5)+ 0.854987973338348); +pos11=find(-2<=T); rigidity(pos11)=10^8*(-0.003748937622488*(T(pos11)+2).^3-0.018449844983174*(T(pos11)+2).^2 -0.057638157095631*(T(pos11)+2)+ 0.746900791092860); + +%Now make sure that rigidity is positive +pos=find(rigidity<0); rigidity(pos)=10^6; Index: /issm/trunk-jpl/src/m/materials/paterson.py =================================================================== --- /issm/trunk-jpl/src/m/materials/paterson.py (revision 13011) +++ /issm/trunk-jpl/src/m/materials/paterson.py (revision 13011) @@ -0,0 +1,50 @@ +from numpy import * + +def paterson(temperature): + + # Local Variables: pos11, pos5, pos10, temperature, pos, T, pos8, pos9, pos6, pos7, pos4, rigidity, pos2, pos3, pos1 + # Function calls: length, zeros, argwhere, paterson, error + #PATERSON - figure out the rigidity of ice for a given temperature + # + # rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). + # temperature is in Kelvin degrees + # + # Usage: + # rigidity=paterson(temperature) + + pos=argwhere(temperature<0.) + if len(pos): + print 'input temperature should be in Kelvin (positive)' + return [] + + T = temperature-273.15 + #The routine below is equivalent to: + # n=3; T=temperature-273; + # %From paterson, + # Temp=[0;-2;-5;-10;-15;-20;-25;-30;-35;-40;-45;-50]; + # A=[6.8*10^-15;2.4*10^-15;1.6*10^-15;4.9*10^-16;2.9*10^-16;1.7*10^-16;9.4* + # 10^-17;5.1*10^-17;2.7*10^-17;1.4*10^-17;7.3*10^-18;3.6*10^-18];;%s-1(kPa-3) + # %Convert into rigidity B + # B=A.^(-1/n)*10^3; %s^(1/3)Pa + # %Now, do a cubic fit between Temp and B: + # fittedmodel=fit(Temp,B,'cubicspline'); + # rigidity=fittedmodel(temperature); + + rigidity=zeros(len(T)) + pos1=argwhere(T<=-45); rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2 -0.325004442485481*(T[pos1]+50)+ 6.524779401948101) + pos2=argwhere(logical_and(-45<=T,T<-40)); rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2 -0.230243014094813*(T[pos2]+45)+ 5.154964909039554) + pos3=argwhere(logical_and(-40<=T,T<-35)); rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+ 0.002886649363879*(T[pos3]+40)**2 -0.179411542205399*(T[pos3]+40)+ 4.149132666831214) + pos4=argwhere(logical_and(-35<=T,T<-30)); rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2 -0.145089762507325*(T[pos4]+35)+ 3.333333333333331) + pos5=argwhere(logical_and(-30<=T,T<-25)); rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2 -0.111773554501713*(T[pos5]+30)+ 2.696559088937191) + pos6=argwhere(logical_and(-25<=T,T<-20)); rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2 -0.088217055680511*(T[pos6]+25)+ 2.199331606342181) + pos7=argwhere(logical_and(-20<=T,T<-15)); rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+ 0.001578771886910*(T[pos7]+20)**2 -0.070194372551690*(T[pos7]+20)+ 1.805165505978111) + pos8=argwhere(logical_and(-15<=T,T<-10)); rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2 -0.044137585824322*(T[pos8]+15)+ 1.510778053489523) + pos9=argwhere(logical_and(-10<=T,T<-5)); rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3- 0.009863871256831*(T[pos9]+10)**2 -0.075294014815659*(T[pos9]+10)+ 1.268434288203714) + pos10=argwhere(logical_and(-5<=T,T<-2)); rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2 -0.048160403003748*(T[pos10]+5)+ 0.854987973338348) + pos11=argwhere(-2<=T); rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2 -0.057638157095631*(T[pos11]+2)+ 0.746900791092860) + + #Now make sure that rigidity is positive + pos=argwhere(rigidity<0); rigidity[pos]=1**6 + + return rigidity + Index: sm/trunk-jpl/src/m/meca/carter.m =================================================================== --- /issm/trunk-jpl/src/m/meca/carter.m (revision 13010) +++ (revision ) @@ -1,28 +1,0 @@ -function attenuation=carter(temperature) -%CARTER - attenuation as a function of temperature -% -% TWO WAY - Attenuation (in dB/m) as a function of temperature (K) -% From Carter at al. 2007 (Radar-based subglacial lake classification in Antarctica) -% Figure 4 -% -% Usage: -% attenuation=carter(temperature) - -if(temperature<0) - error('input temperature should be in Kelvin (positive)'); -end -T=temperature-273.15; - -Temp=[-50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0]'; -A=[0.0015 0.002 0.003 0.0042 0.0055 0.0083 0.012 0.0175 0.026 0.038 0.055]'; - -%Now, do a cubic fit between Temp and B: -[cfun,gof,output]=fit(Temp,A,'cubicspline'); -%breaks=cfun.p.breaks; -%coeff=cfun.p.coefs; - -%Calculate attenuation -attenuation=cfun(T); - -%Make it a 2 way attenuation -attenuation=2*attenuation; Index: /issm/trunk-jpl/src/m/meca/cfl_step.m =================================================================== --- /issm/trunk-jpl/src/m/meca/cfl_step.m (revision 13011) +++ /issm/trunk-jpl/src/m/meca/cfl_step.m (revision 13011) @@ -0,0 +1,23 @@ +function maxtime=cfl_step(md,vx,vy); +%CFL_STEP - return the maximum time step for the model in years +% +% Dt < 0.5 / ( u/Dx +v/Dy ) +% +% Usage: +% maxtime=cfl_step(md,vx,vy); +% +% Example: +% dt=cfl_step(md,md.results.DiagnosticSolution.Vx,md.results.DiagnosticSolution.Vy) + +%Check length of velocities +if size(vx,1)~=md.mesh.numberofvertices & size(vy,1)~=md.mesh.numberofvertices, + error('timestpes error message: size of velocity components must be the same as md.mesh.numberofvertices'); +end + +index=md.mesh.elements; +edgex=max(md.mesh.x(index),[],2)-min(md.mesh.x(index),[],2); +edgey=max(md.mesh.y(index),[],2)-min(md.mesh.y(index),[],2); +vx=max(abs(vx(index)),[],2); +vy=max(abs(vy(index)),[],2); + +maxtime=1/2*min(1./(vx./edgex+vy./edgey)); Index: sm/trunk-jpl/src/m/meca/lithostaticpressure.m =================================================================== --- /issm/trunk-jpl/src/m/meca/lithostaticpressure.m (revision 13010) +++ (revision ) @@ -1,7 +1,0 @@ -function pressure=lithostaticpressure(rho_ice,g,s,z); -%LITHOSTATICPRESSURE - compute the lithostatic pressure -% -% Usage: -% pressure=lithostacticpressure(rho_ice,g,surface,z) - -pressure=rho_ice*g*(s-z); Index: sm/trunk-jpl/src/m/meca/paterson.m =================================================================== --- /issm/trunk-jpl/src/m/meca/paterson.m (revision 13010) +++ (revision ) @@ -1,42 +1,0 @@ -function rigidity=paterson(temperature) -%PATERSON - figure out the rigidity of ice for a given temperature -% -% rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). -% temperature is in Kelvin degrees -% -% Usage: -% rigidity=paterson(temperature) - -if(temperature<0) - error('input temperature should be in Kelvin (positive)'); -end -T=temperature-273.15; - -%The routine below is equivalent to: - -% n=3; T=temperature-273; -% %From paterson, -% Temp=[0;-2;-5;-10;-15;-20;-25;-30;-35;-40;-45;-50]; -% A=[6.8*10^-15;2.4*10^-15;1.6*10^-15;4.9*10^-16;2.9*10^-16;1.7*10^-16;9.4* -% 10^-17;5.1*10^-17;2.7*10^-17;1.4*10^-17;7.3*10^-18;3.6*10^-18];;%s-1(kPa-3) -% %Convert into rigidity B -% B=A.^(-1/n)*10^3; %s^(1/3)Pa -% %Now, do a cubic fit between Temp and B: -% fittedmodel=fit(Temp,B,'cubicspline'); -% rigidity=fittedmodel(temperature); - -rigidity=zeros(length(T),1); -pos1=find(T<=-45); rigidity(pos1)=10^8*(-0.000292866376675*(T(pos1)+50).^3+ 0.011672640664130*(T(pos1)+50).^2 -0.325004442485481*(T(pos1)+50)+ 6.524779401948101); -pos2=find(-45<=T & T<-40); rigidity(pos2)=10^8*(-0.000292866376675*(T(pos2)+45).^3+ 0.007279645014004*(T(pos2)+45).^2 -0.230243014094813*(T(pos2)+45)+ 5.154964909039554); -pos3=find(-40<=T & T<-35); rigidity(pos3)=10^8*(0.000072737147457*(T(pos3)+40).^3+ 0.002886649363879*(T(pos3)+40).^2 -0.179411542205399*(T(pos3)+40)+ 4.149132666831214); -pos4=find(-35<=T & T<-30); rigidity(pos4)=10^8*(-0.000086144770023*(T(pos4)+35).^3+ 0.003977706575736*(T(pos4)+35).^2 -0.145089762507325*(T(pos4)+35)+ 3.333333333333331); -pos5=find(-30<=T & T<-25); rigidity(pos5)=10^8*(-0.000043984685769*(T(pos5)+30).^3+ 0.002685535025386*(T(pos5)+30).^2 -0.111773554501713*(T(pos5)+30)+ 2.696559088937191); -pos6=find(-25<=T & T<-20); rigidity(pos6)=10^8*(-0.000029799523463*(T(pos6)+25).^3+ 0.002025764738854*(T(pos6)+25).^2 -0.088217055680511*(T(pos6)+25)+ 2.199331606342181); -pos7=find(-20<=T & T<-15); rigidity(pos7)=10^8*(0.000136920904777*(T(pos7)+20).^3+ 0.001578771886910*(T(pos7)+20).^2 -0.070194372551690*(T(pos7)+20)+ 1.805165505978111); -pos8=find(-15<=T & T<-10); rigidity(pos8)=10^8*(-0.000899763781026*(T(pos8)+15).^3+ 0.003632585458564*(T(pos8)+15).^2 -0.044137585824322*(T(pos8)+15)+ 1.510778053489523); -pos9=find(-10<=T & T<-5); rigidity(pos9)=10^8*(0.001676964325070*(T(pos9)+10).^3- 0.009863871256831*(T(pos9)+10).^2 -0.075294014815659*(T(pos9)+10)+ 1.268434288203714); -pos10=find(-5<=T & T<-2); rigidity(pos10)=10^8*(-0.003748937622487*(T(pos10)+5).^3+0.015290593619213*(T(pos10)+5).^2 -0.048160403003748*(T(pos10)+5)+ 0.854987973338348); -pos11=find(-2<=T); rigidity(pos11)=10^8*(-0.003748937622488*(T(pos11)+2).^3-0.018449844983174*(T(pos11)+2).^2 -0.057638157095631*(T(pos11)+2)+ 0.746900791092860); - -%Now make sure that rigidity is positive -pos=find(rigidity<0); rigidity(pos)=10^6; Index: sm/trunk-jpl/src/m/meca/paterson.py =================================================================== --- /issm/trunk-jpl/src/m/meca/paterson.py (revision 13010) +++ (revision ) @@ -1,50 +1,0 @@ -from numpy import * - -def paterson(temperature): - - # Local Variables: pos11, pos5, pos10, temperature, pos, T, pos8, pos9, pos6, pos7, pos4, rigidity, pos2, pos3, pos1 - # Function calls: length, zeros, argwhere, paterson, error - #PATERSON - figure out the rigidity of ice for a given temperature - # - # rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). - # temperature is in Kelvin degrees - # - # Usage: - # rigidity=paterson(temperature) - - pos=argwhere(temperature<0.) - if len(pos): - print 'input temperature should be in Kelvin (positive)' - return [] - - T = temperature-273.15 - #The routine below is equivalent to: - # n=3; T=temperature-273; - # %From paterson, - # Temp=[0;-2;-5;-10;-15;-20;-25;-30;-35;-40;-45;-50]; - # A=[6.8*10^-15;2.4*10^-15;1.6*10^-15;4.9*10^-16;2.9*10^-16;1.7*10^-16;9.4* - # 10^-17;5.1*10^-17;2.7*10^-17;1.4*10^-17;7.3*10^-18;3.6*10^-18];;%s-1(kPa-3) - # %Convert into rigidity B - # B=A.^(-1/n)*10^3; %s^(1/3)Pa - # %Now, do a cubic fit between Temp and B: - # fittedmodel=fit(Temp,B,'cubicspline'); - # rigidity=fittedmodel(temperature); - - rigidity=zeros(len(T)) - pos1=argwhere(T<=-45); rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2 -0.325004442485481*(T[pos1]+50)+ 6.524779401948101) - pos2=argwhere(logical_and(-45<=T,T<-40)); rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2 -0.230243014094813*(T[pos2]+45)+ 5.154964909039554) - pos3=argwhere(logical_and(-40<=T,T<-35)); rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+ 0.002886649363879*(T[pos3]+40)**2 -0.179411542205399*(T[pos3]+40)+ 4.149132666831214) - pos4=argwhere(logical_and(-35<=T,T<-30)); rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2 -0.145089762507325*(T[pos4]+35)+ 3.333333333333331) - pos5=argwhere(logical_and(-30<=T,T<-25)); rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2 -0.111773554501713*(T[pos5]+30)+ 2.696559088937191) - pos6=argwhere(logical_and(-25<=T,T<-20)); rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2 -0.088217055680511*(T[pos6]+25)+ 2.199331606342181) - pos7=argwhere(logical_and(-20<=T,T<-15)); rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+ 0.001578771886910*(T[pos7]+20)**2 -0.070194372551690*(T[pos7]+20)+ 1.805165505978111) - pos8=argwhere(logical_and(-15<=T,T<-10)); rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2 -0.044137585824322*(T[pos8]+15)+ 1.510778053489523) - pos9=argwhere(logical_and(-10<=T,T<-5)); rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3- 0.009863871256831*(T[pos9]+10)**2 -0.075294014815659*(T[pos9]+10)+ 1.268434288203714) - pos10=argwhere(logical_and(-5<=T,T<-2)); rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2 -0.048160403003748*(T[pos10]+5)+ 0.854987973338348) - pos11=argwhere(-2<=T); rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2 -0.057638157095631*(T[pos11]+2)+ 0.746900791092860) - - #Now make sure that rigidity is positive - pos=argwhere(rigidity<0); rigidity[pos]=1**6 - - return rigidity - Index: /issm/trunk-jpl/src/m/miscellaneous/isnans.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/isnans.m (revision 13011) +++ /issm/trunk-jpl/src/m/miscellaneous/isnans.m (revision 13011) @@ -0,0 +1,15 @@ +function returnvalue=isnans(array) +%ISNANS: figure out if an array is nan. wrapper to isnan from matlab which stupidly does not allow this test for structures! +% +% Usage: isnans(array) +% +% See also : ISNAN + + +if isstruct(array), + returnvalue=0; +elseif iscell(array) + returnvalue=0; +else + returnvalue=isnan(array); +end Index: /issm/trunk-jpl/src/m/miscellaneous/isnans.py =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/isnans.py (revision 13011) +++ /issm/trunk-jpl/src/m/miscellaneous/isnans.py (revision 13011) @@ -0,0 +1,18 @@ +import numpy + +def isnans(array): + """ + ISNANS: figure out if an array is nan. wrapper to isnan from matlab which stupidly does not allow this test for structures! + + Usage: isnans(array) + + See also : ISNAN + """ + + if isinstance(array,(tuple,list,dict)): + returnvalue=0 + else: + returnvalue=numpy.isnan(array) + + return returnvalue + Index: /issm/trunk-jpl/src/m/miscellaneous/round_ice.m =================================================================== --- /issm/trunk-jpl/src/m/miscellaneous/round_ice.m (revision 13011) +++ /issm/trunk-jpl/src/m/miscellaneous/round_ice.m (revision 13011) @@ -0,0 +1,36 @@ +function new_x=round_ice(x,numnonzeros) +%ROUND_ICE - rounds up x so that it has only numnonzeros non zero digits +% +% numnonzeros must be an integer larger or equal to 1 +% +% Usage: +% new_x=round_ice(x,numnonzeros) + +%some checks +if (nargin ~=2 | nargout>1), + error('round_ice usage: new_x=round_ice(x,numonzeros)'); +end +if ~isnumeric(x) + error('round_ice error message: x must be a number and numzeros an integer'); +end +if round(numnonzeros)~=numnonzeros + error('round_ice error message: numnonzeros must be an integer larger or equal to 1') +end +if any(numnonzeros<1) + error('round_ice error message: numnonzeros must be an integer larger or equal to 1') +end +if (length(numnonzeros)~=1 & size(numnonzeros)~=size(x)) + error('round_ice error message: numnonzeros must be an integer larger or equal to 1 or a list of integers of length length(x)') +end + +%figure out how long x is +lengthx=ceil(log10(abs(x))); + +%if x contains 0, lengthx=-Inf +lengthx(isinf(lengthx))=1; + +%get its sign +si=sign(x); + +%rule out zeros +new_x=si.*round(abs(x).*10.^(-lengthx+numnonzeros)).*10.^(lengthx-numnonzeros);