pro calcium

;calcium derivation program, for experimental records

;use calciumsim for simulations

dummy=''

contout=''

outnamo=''

fnama=''

;IMPORTANT PARAMETERS

zoom=3.5 & laser=3.& gain=1000.

kmax = 2;1.45

ca0=0.05

;DEFINITIONS

nwi=44 & nwm=nwi-1 & nmm=nwi-2 & nmmm=nwi-3 & nwp=nwi+1

nl = 768

;dimension definitions ********************

jump0:

nwih=fix(nwi/2)

caimg=fltarr(nwi,nwp)

bb=fltarr(nwi)

dyarr=fltarr(nwp,nwp) & dyrr=fltarr(nwi,nwp)

z=fltarr(nwp,nwp) & ri=fltarr(nwi,nwp)

dzdt=fltarr(nwi,nwp) & dzdx=fltarr(nwi,nwp) & d2z=fltarr(nwi,nwp)

avgarr=fltarr(nwi,nwp)

flav=fltarr(nwp,nwp) & dyav=fltarr(nwp) & dccd=fltarr(nwi,nwp)



; UNITS ARE MICROMOLAR, MILLISECOND, MICRON, PICOAMPERES

dcd= 0.02 ;diff. coef. dye. (in micron2/ms)

dcc=0.35 ;diff. coef. Ca

;rate constants are in per micromolar per millisecond (multiply the usual /M/s by 10^-9)

fkd =1 & fkon=0.032;1.03 & fkon=0.032

fkof=fkon*fkd ;off rate const fluo

;miscellaneous

dt=2.0

dx= 0.500/zoom

r = 100.0 ;Fmax/Fmin value for lower estimate of Ca

;INTERPOLATION to go from an array of 45 pixels in width to one of 44.

;Additional goal is not to have x=0 (singular point)

rin=float(nwi)/float(nwm) & dx=dx*rin ;used in interpolation

x=findgen(nwp)

x=interpolate(x,rin*findgen(nwi))

x=x-nwi/2

x=rebin(x,nwi,nwp)*dx

;& xd=dx*findgen(44) ;makes the interpolated x coordinate

;td=dt*findgen(44) ;makes the t coordinate

dyav=0 ;used for averaging dye

;IMAGE INPUT, ENTERS AN AVERAGE SPARK, MADE BY ANOTHER PROGRAM. AN ARRAY OF 45 PIXEL WIDE BY 46 LONG, [DYE] IN 1ST ROW.

read, fnama, prompt='Enter average file,( no ext., default .avg): '

fname=fnama+'.avg'

fnama=STRMID(fnama,0,7);NOTE REDEFINITION FOR SHORTER FNAMA

jump12:

counti=1

openu,u,fname,/get_lun & zu=assoc(u,fltarr(nwp,nwp+1))

avgarr=zu(counti-1) ;average spark, with dye data in first row

close,u & free_lun,u

gain=avgarr(nwi,nwp)*1000. ;gain is recovered, was stored in average spark

effgain=exp((gain-1000.)*0.011) ;ratio of actual gain and reference gain

;display of average

ah=-0.75 & bh=0.25 & ch=1.0 & back=255

loadct,0

set_shading,light=[ah,bh,ch]

anglex=45

anglez=45

jump07:

shade_surf,avgarr(*,1:45),ax=anglex,az=anglez,charsize=2,xstyle=4,ystyle =4, background = back, color =1 ,zstyle=1;zrange=[0,200],

read,dummy,prompt='any key, b to change background, angles:'

if dummy ne 'b' then goto,jump08

read,back,anglez,anglex,prompt='enter background level, anglez, anglex:'

goto,jump07

jump08:

;***************

flav=avgarr(*,1:nwp) ;average spark fluorescence, omitting 1st row

dyav=avgarr(*,0) ;dye, a function of x, called dyav 'cause this is always an average over many sparks

fluo=total(dyav)/n_elements(dyav) ;the average of dye over x

;will construct dyarr, an array of total dye (x) , a constant in time, symmetric around spark center

for j=0,nwi do dyav(j)=(dyav(j)+dyav(nwi-j))/2 ;symmetrization

for j=0,nwi do dyarr(j,*)=dyav(j) ;converted to a 45 x 45 array, constant in time

; will regenerate baseline fluorescence, a function of x

bb =laser* dyav*kmax*exp((gain-1000.)*0.011)*ca0/fkd

bbo=total(bb)/n_elements(bb) ;average baseline fluorescence

;bbl=rebin(flav(0:nwi,0:9),nwp,1) ;local, used in other programs



fluo=total(dyav)/n_elements(dyav) ;average dye concentration over baseline

dcd_x=dcd*total(bb)/bb/n_elements(bb) ;dye diffusion constant corrected for a local dye binding factor,

; which explains unequal dye distribution; in re9 it's done with bbl

dcd_x=rebin(dcd_x,45,45) ;made to an array

;important calculations start here, with Ca-bound dye



for i=0,nwi do begin ;i is boxcar time

for j=0,nwi do begin ;j varies over x coordinate

z(j,i)=(r*(flav(j,i))/effgain/laser-dyav(j)*kmax)/(r-1)/kmax

;z is Ca-bound dye, calculated straightforwardly from fluorescence

endfor

endfor

;***********************

for i=0,nwi do begin ;loop over time

ri(*,i)=interpolate(z(*,i),rin*findgen(nwi)) ;converts to 44 pixels instead of 45

dyrr(*,i)=interpolate(dyarr(*,i),rin*findgen(nwi)) ;same for dye array

dccd(*,i)=interpolate(dcd_x(*,i),rin*findgen(nwi)) ;same for dye binding constant array

endfor

z=ri ;z is [dye:Ca](x,t) after interpolation

dyarr=dyrr ;dyarr is [dye total](x,t), a consant in t but not in x

dertb,z,dzdt,dt ;subroutine calculates dzdt, time derivative of Ca-bound dye

dzdt(*,[0,1,nmm,nwm])=0.0 ;making edges nice, probably inadequate for large sparks

dzdt([0,1,nmm,nwm],*)=0.0

derxb,z,dzdx,dx ;subroutine calculates dzdx, first space derivative

dzdx(*,[0,1,nmm,nwm])=0.0 & dzdx([0,1,nmm,nwm],*)=0.0

;calculate second derv with space

derxb,dzdx,d2z,dx ;2nd derivative in d2z

d2z(*,[0,1,nmm,nwm])=0.0 & d2z([0,1,nmm,nwm],*)=0.0

diff_term=-dccd*(d2z+2*dzdx/x) ;laplacian term in [Ca2+] calculation, assumed spherically symmetric

dye_term=diff_term+dzdt ;note all dye-related terms added

;*************************

caimg=((dye_term+fkof*z)/fkon/(dyarr-z))>ca0 ;dynamically calculated [Ca2+]

for i=0,nwm do begin

for j=0,nwi do begin

caimg(i,j)=(caimg(i,j)+caimg(nwm-i,j))/2 ;symmetrizing, just in case

endfor

endfor

;DISPLAY AND OUTPUT *************************

result=congrid(caimg,360,360,cubic=-0.5)

loadct,4

tvscl,result

read,dummy,prompt='any key, s to stop'

if dummy eq 's' then stop

loadct,0

set_shading,light=[ah,bh,ch]

shade_surf,caimg(1:nwi-2,10:nwi-2),ax=45,az=45,color=1, xstyle=4,ystyle=4,charsize=2,background = 255

read, dummy, prompt='f for new file, any key to stop:'

if dummy eq 'f' then goto, jump0

stop

end