/* * Copyright (C) 2013-2015, Gaetan Bisson . * * Permission to use, copy, modify, and/or distribute this software for any * purpose with or without fee is hereby granted, provided that the above * copyright notice and this permission notice appear in all copies. * * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY * SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION * OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ /* * Decox. Decreasing continuous exponentials. * * Decox merely aims to give its author a better understanding of the behavior * of decreasing exponentials in the context of Haldanian decompression models. * This software is unfit for actual dive planning; never rely on its output! * * Units are bars, liters, meters and minutes. * * The "gas" argument to many procedures is a function that takes as input a * depth and returns the partial pressure of nitrogen the diver is exposed to. * For planning, "mix" triplets consist of a label for, the maximum operating * depth of, and the gas itself. Helper functions are provided for open circuit * mixes and closed circuit setpoints. * * Run with PARI/GP: * * gp -q decox.gp */ /* **************************************************************************** * * GLOBAL PHYSICAL CONSTANTS * */ /* infinity for physicists */ MAX=9999; /* atmosphere pressure */ atm=1.013; /* saltwater density */ density=1.03; /* time step for discretization */ step=0.01; /* compartments' half-time and m-value from ZH-L16B */ ht=[5, 8, 12.5, 18.5, 27, 38.3, 54.3, 77, 109, 146, 187, 239, 305, 390, 498, 635]; m0=[29.6, 25.4, 22.5, 20.3, 19.0, 17.5, 16.5, 15.7, 15.2, 14.6, 14.2, 13.9, 13.4, 13.2, 12.9, 12.7]; md=[1.7928, 1.5352, 1.3847, 1.2780, 1.2306, 1.1857, 1.1504, 1.1223, 1.0999, 1.0844, 1.0731, 1.0635, 1.0552, 1.0478, 1.0414, 1.0359]; /* gradient factor coefficients and pressure points */ gf1c=0.9; gf2c=0.3; gf1p=1; gf2p=7; /* safe ascent rate */ rate=9; /* surface air consumption */ sac=15; /* nitrogen partial pressure at depth */ ccr(sp)=(depth)->max(0,pressure(depth)-sp); ean(o2)=(depth)->(1-o2/100)*pressure(depth); air=ean(21); /* mix triplets further include a label and MOD */ CCR(sp)=[Strprintf("SP=%4.2f",sp),MAX,ccr(sp)]; EAN(o2)=[Strprintf("EAN%i",o2),(1600/o2-10)\3*3,ean(o2)]; AIR=EAN(21);AIR[1]="AIR"; /* **************************************************************************** * * DECOMPRESSION FUNCTIONS AND CONSTANTS * * Only two global variables are used for planning: * - run, the current runtime of the plan * - V, a map of the current gas volumes */ /* compartments' indices */ N=[1..#ht]; /* ambient pressure at depth */ pressure(depth)=atm+density*depth/10; /* gradient factor at ambiant ressure */ gf(p)=(gf1c*(gf2p-p)-gf2c*(gf1p-p))/(gf2p-gf1p); /* m-value from ambiant pressure */ mval=[(p)->p+gf(p)*(m0[i]/10+md[i]*(p-1)-p) | i<-N]; /* nitrogen tension vector after time breathing pn */ update(pp,time,pn)=[pn+0.5^(time/ht[i])*(pp[i]-pn) | i<-N]; /* time breathing pn until any compartment greater than lim */ time2any(pp,pn,lim)={ my(nhl=(pn,pp,lim)->q=(pn-pp)/(pn-lim);if(q>0,log(q)/log(2),if(pn>lim,0,MAX))); vecmin([ht[i]*nhl(pn,pp[i],lim[i]) | i<-N]); } /* time breathing pn until all compartments smaller than lim */ time2all(pp,pn,lim)={ my(nhl=(pn,pp,lim)->q=(pn-pp)/(pn-lim);if(q>0,log(q)/log(2),if(pn>lim,MAX,0))); vecmax([ht[i]*nhl(pn,pp[i],lim[i]) | i<-N]); } /* time at surface before 0.7 atm allowable */ nofly(pp)=time2all(pp,air(0),[mval[i](0.7*atm) | i<-N]); /* time at surface before complete desaturation */ desat(pp)=time2all(pp,air(0),[air(0)*1.01 | i<-N]); /* time at depth "from" until direct ascent to "to" allowable */ time2asc(pp,from,to,gas)={ my(depth,worst,stay,stop); if(from0,return(-stay)); if(stop>0,return(stop)); return(-MAX); } /* no decompression limit breathing gas at depth */ ndl(depth,gas)=max(0,-time2asc([air(0) | i<-N],depth,0,gas)); /* mix volume helper functions */ volini()=V=Map(); voladd(s,v)={my(c=0);mapisdefined(V,s,&c);mapput(V,s,c+v);} mixout(s)=if(Vec(s)[1]!="S",my(v=mapget(V,s));printf(" used %4iL = %3ib of 12L %s\n",ceil(v),ceil(v/12),s)); volout()=apply(mixout,Vec(V)); /* decompression schedule segment from depth "from" to "to" */ segment(pp,from,to,mix)={ my(depth,str,gas,up3m,time); depth=from; str=mix[1]; gas=mix[3]; while(depth>to, up3m=max((depth-1)\3*3,to); time=ceil(time2asc(pp,depth,up3m,gas)); if(time>0, pp=update(pp,time,gas(depth)); voladd(str,time*sac*pressure(depth)); run+=time; printf("Stop at %2im for %2i' (%3i)\n",depth,time,run); ); while(depth>up3m, pp=update(pp,step,gas(depth)); voladd(str,step*sac*pressure(depth)); depth-=step*rate; run+=step; ); depth=up3m; ); pp; } /* decompression schedule after [depth,time,mix] vector using mixes from M */ plan(DTM,M)={ my(pp,depth,time,mix,str,gas,tts,mod); M=concat(M,[[0,0,0]]); pp=[air(0) | i<-N]; volini(); run=0; for(i=1,#DTM, depth=DTM[i][1]; time=DTM[i][2]; mix=DTM[i][3]; str=mix[1]; gas=mix[3]; run+=time; pp=update(pp,time,gas(depth)); voladd(str,time*sac*pressure(depth)); printf("Stay at %2im for %2i' (%3i) %s\n",depth,time,run,DTM[i][3][1]); if(i<#DTM,pp=segment(pp,depth,DTM[i+1][1],mix)); ); tts=run; for(j=1,#M, mod=M[j][2]; if (mod>=depth,next); mix=M[j-1]; str=mix[1]; gas=mix[3]; if(j==1,printf("==> No bailout available.\n\n");return); printf("Up from %2im (%3i) %s\n",depth,run,M[j-1][1]); pp=segment(pp,depth,mod,mix); depth=mod; ); printf(" total ascent: %2i'\n",run-tts); volout(); print(); } /* decompression schedule from any point during the dive */ aplan(DTM,M)={ for(i=1,#DTM, dtm=vector(i,j,DTM[j]); plan(dtm,M); ); } /* **************************************************************************** * * USER "INTERFACE" * */ { printf("Decox --- Decreasing continuous exponentials --- Version 7. Copyright (C) 2013-2015, Gaetan Bisson. All rights reserved. See the source code for licensing terms and documentation. Example: mix=CCR(1.2); bottom=[[40,50,mix],[20,50,mix]]; plan(bottom,[mix]); aplan(bottom,[AIR,EAN(50)]); "); } default(timer,0); default(logfile,"~/.decox.log"); default(log,1);