From 71f8dd0dca990b8c1d4bb01f3dae68e1e90f771e Mon Sep 17 00:00:00 2001 From: TimothyADavis Date: Wed, 28 Oct 2015 12:23:55 +0000 Subject: [PATCH 1/5] Function for computing INCLOUDS --- KinMS.pro | 69 ++++++++++++++++++++++++++----------------------------- 1 file changed, 33 insertions(+), 36 deletions(-) diff --git a/KinMS.pro b/KinMS.pro index d605c31..1f316ac 100644 --- a/KinMS.pro +++ b/KinMS.pro @@ -176,7 +176,26 @@ ; ;############################################################################## -pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=sbrad,velrad=velrad,velprof=velprof,filename=galname,diskthick=diskthick,cleanout=cleanout,ra=ra,dec=dec,nsamps=nsamps,cubeout=cubeout,posang=posang,intflux=intflux,inclouds=inclouds,vlos_clouds=vlos_clouds,flux_clouds=flux_clouds,vsys=vsys,restfreq=restfreq,phasecen=phasecen,voffset=voffset,fixseed=fixseed,vradial=vradial,vphasecen=vphasecen,vposang=vposang + +function kinms_samplefromarbdist_onesided,sbrad,sbprof,nsamps,seed,r_flat=r_flat,diskthick=diskthick + px=dblarr(n_elements(sbprof)) + sbprofboost=sbprof*(2*!dpi*abs(sbrad)) ;; boost flux in outer parts to compensate for spreading over 2pi radians + for i=1,n_elements(sbrad)-1 do px[i] = INT_TABULATED((sbrad[0:i]),(sbprofboost[0:i]),/double) ;; integrate surf brightness w.r.t radius + px/=max(px) ;; normalize + pick=randomu(seed[0],nsamps) ;; pick random y to transform + phi=randomu(seed[1],nsamps)*2*!dpi ;; random angle + r_flat=interpol(sbrad,px,pick) ;; invert to get correct distribution + if n_elements(diskthick) gt 1 then diskthick_here=interpol(diskthick,sbrad,r_flat) else diskthick_here=diskthick + zpos=(diskthick_here)*((randomu(seed[2],nsamps)*2)-1) ;; do disk thickness + r_3d = sqrt((r_flat^2)+zpos^2) ;; 3d distance to each + theta=acos(zpos/r_3d) ;; angle out of plane + xpos=((r_3d*cos(phi)*sin(theta))) ;; go to cartesian + ypos=((r_3d*sin(phi)*sin(theta))) ;; go to cartesian + return,TRANSPOSE([[xpos],[ypos],[zpos]]) ;; return INCLOUDS format +end + + +pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=sbrad,velrad=velrad,velprof=velprof,filename=galname,diskthick=diskthick,cleanout=cleanout,ra=ra,dec=dec,nsamps=nsamps,cubeout=cubeout,posang=posang,intflux=intflux,inclouds=inclouds,vlos_clouds=vlos_clouds,flux_clouds=flux_clouds,vsys=vsys,restfreq=restfreq,phasecen=phasecen,voffset=voffset,fixseed=fixseed,vradial=vradial,vphasecen=vphasecen,vposang=vposang,sb_sampfunc=sb_sampfunc ;!EXCEPT = 2 @@ -212,7 +231,7 @@ pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=s if keyword_set(fixseed) then begin seed=indgen(4)+100 ;; fix seeds endif else seed=randomu(seedit,4)*100. ;; random seeds - + if not keyword_set(sb_sampfunc) then sb_sampfunc="kinms_samplefromarbdist_onesided" ;;;; @@ -225,37 +244,17 @@ pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=s vphasecent=(vphasecen-phasecen)/[dx,dy] ;;;; + if not keyword_set(inclouds) then inclouds=call_FUNCTION(sb_sampfunc,sbrad,sbprof,nsamps,seed,r_flat=r_flat,diskthick=diskthick) ;; use the function in sb_sampfunc to setup inclouds, if its not already set -;;;; if using INCLOUDS then set this up - if keyword_set(inclouds) then begin - xpos=reform(INCLOUDS[0,*]/dx) - ypos=reform(INCLOUDS[1,*]/dy) - zpos=reform(INCLOUDS[2,*]/dx) - r_flat=sqrt(((xpos-(phasecen[0]/dx))^2) + ((ypos-(phasecen[1]/dy))^2)) - if keyword_set(VLOS_CLOUDS) then los_vel=VLOS_CLOUDS - endif - + +;;;; set up INCLOUDS for use + xpos=reform(INCLOUDS[0,*]/dx) + ypos=reform(INCLOUDS[1,*]/dy) + zpos=reform(INCLOUDS[2,*]/dx) + if n_elements(r_flat) eq 0 then r_flat=sqrt(((xpos-(phasecen[0]/dx))^2) + ((ypos-(phasecen[1]/dy))^2)) + if keyword_set(VLOS_CLOUDS) then los_vel=VLOS_CLOUDS + -;;;; Else draw from arb. surface brightness distribution in SBPROF ;;; - if not keyword_set(inclouds) then begin - sbprofs=sbprof ;; save originals for later - sbrads=sbrad ;;save originals for later - px=dblarr(n_elements(sbprof)) - sbrad=temporary(sbrad)/dx - sbprof=temporary(sbprof)*(2*!dpi*abs(sbrad)) ;; boost flux in outer parts to compensate for spreading over 2pi radians - for i=1,n_elements(sbrad)-1 do px[i] = INT_TABULATED((sbrad[0:i]),(sbprof[0:i]),/double) ;; integrate surf brightness w.r.t radius - px/=max(px) ;; normalize - pick=randomu(seed[0],nsamps) ;; pick random y to transform - phi=randomu(seed[1],nsamps)*2*!dpi ;; random angle - r_flat=interpol(sbrad,px,pick) ;; invert to get correct distribution - if n_elements(diskthick) gt 1 then diskthick_here=interpol(diskthick,sbrad,r_flat) else diskthick_here=diskthick - zpos=(diskthick_here/dx)*((randomu(seed[2],nsamps)*2)-1) ;; do disk thickness - r_3d = sqrt((r_flat^2)+zpos^2) ;; 3d distance to each - theta=acos(zpos/r_3d) ;; angle out of plane - xpos=((r_3d*cos(phi)*sin(theta))) ;; go to cartesian - ypos=((r_3d*sin(phi)*sin(theta))) ;; go to cartesian - ;;;; - endif ;;;; create velocity structure ;;;; @@ -317,8 +316,9 @@ if not keyword_set(inclouds) then begin INCLOUDS[0,*]=xpos*dx INCLOUDS[1,*]=ypos*dy INCLOUDS[2,*]=zpos*dx - vlos_clouds=los_vel endif +if not keyword_set(vlos_clouds) then vlos_clouds=los_vel + ;;;; ; stop @@ -373,10 +373,7 @@ endif cube*=((INT_TABULATED( sbrad, sbprof))/((total(cube)*dv)/(sqrt(beamsize[0]*beamsize[1])^2))) ;; normalize to get same flux as input sb profile variable endelse endelse -if keyword_set(sbprof) then begin - sbprof=sbprofs ;; restore originals - sbrad=sbrads ;;restore originals -endif + ;;;; write cube ;;;; if keyword_set(galname) then begin From 04f7aa149508efcecfd29b2821aef7d1202d4d30 Mon Sep 17 00:00:00 2001 From: TimothyADavis Date: Wed, 28 Oct 2015 13:39:11 +0000 Subject: [PATCH 2/5] Function for computing velfield --- KinMS.pro | 100 +++++++++++++++++++++++++++++------------------------- 1 file changed, 54 insertions(+), 46 deletions(-) diff --git a/KinMS.pro b/KinMS.pro index 1f316ac..cd486a9 100644 --- a/KinMS.pro +++ b/KinMS.pro @@ -178,6 +178,10 @@ function kinms_samplefromarbdist_onesided,sbrad,sbprof,nsamps,seed,r_flat=r_flat,diskthick=diskthick + + ;;;; This function takes the input radial distribution and draws + ;;;; n samples from under it. It also accounts for disk thickness + px=dblarr(n_elements(sbprof)) sbprofboost=sbprof*(2*!dpi*abs(sbrad)) ;; boost flux in outer parts to compensate for spreading over 2pi radians for i=1,n_elements(sbrad)-1 do px[i] = INT_TABULATED((sbrad[0:i]),(sbprofboost[0:i]),/double) ;; integrate surf brightness w.r.t radius @@ -194,8 +198,45 @@ function kinms_samplefromarbdist_onesided,sbrad,sbprof,nsamps,seed,r_flat=r_flat return,TRANSPOSE([[xpos],[ypos],[zpos]]) ;; return INCLOUDS format end +function kinms_create_velfield_onesided,velrad,velprof,r_flat,inc,posang,gassigma,seed,xpos,ypos,vphasecent,vposang=vposang,vradial=vradial,inc_rad=inc_rad,posang_rad=posang_rad + + + ;;;; This function takes the input circular velocity distribution + ;;;; and creates the velocity field taking into account warps, + ;;;; inflow/outflow etc as required by the keywords passed. + + vrad=interpol(velprof,velrad,r_flat) + los_vel=dblarr(n_elements(vrad)) + + ;;; in case of warps/inflow/outflow then setup vectors to deal with this + if n_elements(vradial) gt 1 then vradial_rad=interpol(vradial,velrad,r_flat) else vradial_rad=fltarr(n_elements(r_flat))+vradial + if n_elements(vposang) gt 1 then vposang_rad=interpol(vposang,velrad,r_flat) else vposang_rad=vposang + ;;;; + + ;;;;; include random motions + if keyword_set(gassigma) then begin + veldisp=randomn(seed[3],n_elements(xpos)) + if n_elements(gassigma) gt 1 then veldisp=temporary(veldisp)*interpol(gassigma,velrad,r_flat) else veldisp=temporary(veldisp)*gassigma + endif else veldisp=fltarr(n_elements(xpos)) + ;;;;; + + ;;;;; project velocities taking into account inclination, which may change with radius + w=where(xpos+vphasecent[0] ne 0.0 and ypos+vphasecent[1] ne 0.0) + ang2rot=(90-(posang_rad-vposang_rad)) + los_vel[w]=veldisp[w] ;; random motions + los_vel[w]+=vrad[w]*(cos(atan(float(ypos[w]+vphasecent[1])/float(xpos[w]+vphasecent[0]))+(!dtor*(90-ang2rot[w])))*sin(inc_rad[w]*!dtor)) ;; circular motions + los_vel[w]+=vradial_rad[w]*(sin(atan(float(ypos[w]+vphasecent[1])/float(xpos[w]+vphasecent[0]))+(!dtor*(90-ang2rot[w])))*sin(inc_rad[w]*!dtor)) ;; inflow/outflow + w=where(xpos+vphasecent[0] eq 0.0 and ypos+vphasecent[1] eq 0.0) + if w[0] ne -1 then los_vel[w]=veldisp[w]+(vrad[w]*sin(inc_rad[w]*!dtor))+ vradial_rad[w]*(sin(atan(float(ypos[w])/float(xpos[w]))+(!dtor*(90-ang2rot[w])))*sin(inc_rad[w]*!dtor)) + w=where(xpos+vphasecent[0] gt 0.0) + los_vel[w]=temporary(los_vel[w])*(-1) + ;;;;;; + + return,los_vel + end + -pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=sbrad,velrad=velrad,velprof=velprof,filename=galname,diskthick=diskthick,cleanout=cleanout,ra=ra,dec=dec,nsamps=nsamps,cubeout=cubeout,posang=posang,intflux=intflux,inclouds=inclouds,vlos_clouds=vlos_clouds,flux_clouds=flux_clouds,vsys=vsys,restfreq=restfreq,phasecen=phasecen,voffset=voffset,fixseed=fixseed,vradial=vradial,vphasecen=vphasecen,vposang=vposang,sb_sampfunc=sb_sampfunc +pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=sbrad,velrad=velrad,velprof=velprof,filename=galname,diskthick=diskthick,cleanout=cleanout,ra=ra,dec=dec,nsamps=nsamps,cubeout=cubeout,posang=posang,intflux=intflux,inclouds=inclouds,vlos_clouds=vlos_clouds,flux_clouds=flux_clouds,vsys=vsys,restfreq=restfreq,phasecen=phasecen,voffset=voffset,fixseed=fixseed,vradial=vradial,vphasecen=vphasecen,vposang=vposang,sb_sampfunc=sb_sampfunc,vel_func=vel_func ;!EXCEPT = 2 @@ -232,7 +273,7 @@ pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=s seed=indgen(4)+100 ;; fix seeds endif else seed=randomu(seedit,4)*100. ;; random seeds if not keyword_set(sb_sampfunc) then sb_sampfunc="kinms_samplefromarbdist_onesided" - + if not keyword_set(vel_func) then vel_func="kinms_create_velfield_onesided" ;;;; ;;;; work out images sizes ;;;; @@ -244,8 +285,10 @@ pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=s vphasecent=(vphasecen-phasecen)/[dx,dy] ;;;; - if not keyword_set(inclouds) then inclouds=call_FUNCTION(sb_sampfunc,sbrad,sbprof,nsamps,seed,r_flat=r_flat,diskthick=diskthick) ;; use the function in sb_sampfunc to setup inclouds, if its not already set - + if not keyword_set(inclouds) then begin + inclouds=call_FUNCTION(sb_sampfunc,sbrad,sbprof,nsamps,seed,r_flat=r_flat,diskthick=diskthick) ;; use the function in sb_sampfunc to setup inclouds, if its not already set + r_flat/=dx + endif ;;;; set up INCLOUDS for use xpos=reform(INCLOUDS[0,*]/dx) @@ -259,44 +302,19 @@ pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=s ;;;; create velocity structure ;;;; if not keyword_set(VLOS_CLOUDS) then begin - velrad=temporary(velrad)/dx - vrad=interpol(velprof,velrad,r_flat) - los_vel=dblarr(n_elements(vrad)) - ;;; in case of warps/inflow/outflow then setup vectors to deal - ;;; with this - if n_elements(inc) gt 1 then inc_rad=interpol(inc,sbrad,r_flat) else inc_rad=fltarr(n_elements(r_flat))+inc - if n_elements(vradial) gt 1 then vradial_rad=interpol(vradial,sbrad,r_flat) else vradial_rad=fltarr(n_elements(r_flat))+vradial - if n_elements(posang) gt 1 then posang_rad=interpol(posang,sbrad,r_flat) else posang_rad=posang - if n_elements(vposang) gt 1 then vposang_rad=interpol(vposang,sbrad,r_flat) else vposang_rad=vposang - ;;;; - ;;;;; include random motions - if keyword_set(gassigma) then begin - veldisp=randomn(seed[3],n_elements(xpos)) - if n_elements(gassigma) gt 1 then veldisp=temporary(veldisp)*interpol(gassigma,sbrad,r_flat) else veldisp=temporary(veldisp)*gassigma - endif else veldisp=fltarr(n_elements(xpos)) - ;;;;; - ;;;;; project velocities taking into account inclination, which may change with radius - w=where(xpos+vphasecent[0] ne 0.0 and ypos+vphasecent[1] ne 0.0) - ang2rot=(90-(posang_rad-vposang_rad)) - los_vel[w]=veldisp[w] ;; random motions - los_vel[w]+=vrad[w]*(cos(atan(float(ypos[w]+vphasecent[1])/float(xpos[w]+vphasecent[0]))+(!dtor*(90-ang2rot[w])))*sin(inc_rad[w]*!dtor)) ;; circular motions - los_vel[w]+=vradial_rad[w]*(sin(atan(float(ypos[w]+vphasecent[1])/float(xpos[w]+vphasecent[0]))+(!dtor*(90-ang2rot[w])))*sin(inc_rad[w]*!dtor)) ;; inflow/outflow - w=where(xpos+vphasecent[0] eq 0.0 and ypos+vphasecent[1] eq 0.0) - if w[0] ne -1 then los_vel[w]=veldisp[w]+(vrad[w]*sin(inc_rad[w]*!dtor))+ vradial_rad[w]*(sin(atan(float(ypos[w])/float(xpos[w]))+(!dtor*(90-ang2rot[w])))*sin(inc_rad[w]*!dtor)) - w=where(xpos+vphasecent[0] gt 0.0) - los_vel[w]=temporary(los_vel[w])*(-1) -;;;; + + if n_elements(inc) gt 1 then inc_rad=interpol(inc,velrad,r_flat) else inc_rad=fltarr(n_elements(r_flat))+inc + if n_elements(posang) gt 1 then posang_rad=interpol(posang,velrad,r_flat) else posang_rad=posang -;;;; project face on clouds to desired inclination ;;;; + los_vel=call_FUNCTION(vel_func,velrad/dx,velprof,r_flat,inc,posang,gassigma,seed,xpos,ypos,vphasecent,vposang=vposang,vradial=vradial,inc_rad=inc_rad,posang_rad=posang_rad) + ;;;; project face on clouds to desired inclination ;;;; c = cos(inc_rad/!radeg) s = sin(inc_rad/!radeg) x2 = xpos y2 = c*ypos + s*zpos - z2 = -s*ypos + c*zpos - + z2 = -s*ypos + c*zpos if n_elements(posang) ge 1 then begin ;; rotate gas on sky to required angle - ang=90-posang_rad c = cos(ang/!radeg) s = sin(ang/!radeg) @@ -310,18 +328,8 @@ pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=s y2=ypos z2=zpos endelse -;;;; return the clouds used incase they are useful in future calls ;;; -if not keyword_set(inclouds) then begin - inclouds=dblarr(3,n_elements(xpos)) - INCLOUDS[0,*]=xpos*dx - INCLOUDS[1,*]=ypos*dy - INCLOUDS[2,*]=zpos*dx -endif -if not keyword_set(vlos_clouds) then vlos_clouds=los_vel + if not keyword_set(vlos_clouds) then vlos_clouds=los_vel -;;;; - - ; stop ;;;; now add the flux into the cube From a7e798b982b3fd75aa44392f20ae660822be1832 Mon Sep 17 00:00:00 2001 From: TimothyADavis Date: Wed, 28 Oct 2015 15:18:45 +0000 Subject: [PATCH 3/5] Updated documentation --- KinMS.pro | 43 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 36 insertions(+), 7 deletions(-) diff --git a/KinMS.pro b/KinMS.pro index cd486a9..ad9f5b2 100644 --- a/KinMS.pro +++ b/KinMS.pro @@ -22,7 +22,7 @@ ; KinMS ; ; PURPOSE: -; Simulate rotating molecular gas distributions, and produce +; Simulate rotating gas distributions, and produce ; datacubes in FITS format. ; ; CALLING SEQUENCE: @@ -46,7 +46,7 @@ ; INC Inclination angle of the gas disc on the sky ; (degrees). Can input a constant or a vector, ; giving the inclination as a function of the -; radius vector SBRAD (in order to model warps etc) +; radius vector VELRAD (in order to model warps etc) ; ; ; @@ -132,20 +132,46 @@ ; FIXSEED Fix the seeds supplied to the random number generators ; so that the output will always be the same for given input ; parameters -; +; VRADIAL Magnitude of inflow/outflowing motions (km/s). Negative +; numbers here are inflow, positive numbers denote +; outflow. These are included in the velocity field using +; formalism of KINEMETRY (Krajnović et al. 2006 MNRAS, 366, 787). +; Can input a constant or a vector, giving the radial +; motion as a function of the radius vector +; VELRAD. Default is no inflow/outflow. +; SB_SAMPFUNC The function to be used to generate a cloud of point +; sources from the input surface brightness +; profile. Default is to draw from a one sided +; probability distribution. This can be swapped out by +; using a user designed function which accepts the +; same inputs/produces the same outputs if one has a +; different behavior in mind. +; VEL_FUNC The function to be used to generate the velocity +; field for the cloud of point sources in INCLOUDS +; based on the input circular velocity profile. +; Default is a radially symetric velocity curve +; This can be swapped out by using a user +; designed function which accepts the same +; inputs/produces the same outputs. + + ; OPTIONAL OUTPUTS: ; ; CUBEOUT Returns the created cube as a 3 dimensional vector ; +; INCLOUDS Returns the cloud of point sources created for your +; model, incase they are useful for a future call. ; +; VLOS_CLOUDS Returns the velocity of each point source created for your +; model, incase they are useful for a future call. +; +; ; SIDE EFFECTS: ; The output file is overwritten if it exists. ; ; RESTRICTIONS: ; Requires IDL 5.0 or higher (square bracket array -; syntax). Requires the latest IDL astrolib, along with the -; hist_nd routine (a version of which is distributed with the -; package for ease). +; syntax). Requires the latest IDL astrolib. ; ; ; USAGE EXAMPLE: @@ -173,6 +199,8 @@ ; of inflow and outflow ; 02/09/2015 v1.5 - Update to handling of position angles, added ; seperate velocity position angle/phasecentre +; 28/10/2015 v1.6 - Updated to allow hot swapping of SBprofile and +; velprofile handling. ; ;############################################################################## @@ -238,7 +266,8 @@ function kinms_create_velfield_onesided,velrad,velprof,r_flat,inc,posang,gassigm pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=sbrad,velrad=velrad,velprof=velprof,filename=galname,diskthick=diskthick,cleanout=cleanout,ra=ra,dec=dec,nsamps=nsamps,cubeout=cubeout,posang=posang,intflux=intflux,inclouds=inclouds,vlos_clouds=vlos_clouds,flux_clouds=flux_clouds,vsys=vsys,restfreq=restfreq,phasecen=phasecen,voffset=voffset,fixseed=fixseed,vradial=vradial,vphasecen=vphasecen,vposang=vposang,sb_sampfunc=sb_sampfunc,vel_func=vel_func ;!EXCEPT = 2 - + + ;;;; Main procedure ;;;; set defaults ;;;; if not keyword_set(nsamps) then nsamps=10000 From f0cbf3e0279e1b6aad653eeef3fa805e5f0003dd Mon Sep 17 00:00:00 2001 From: TimothyADavis Date: Wed, 28 Oct 2015 15:29:22 +0000 Subject: [PATCH 4/5] Changed inclouds format for speed --- KinMS.pro | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/KinMS.pro b/KinMS.pro index ad9f5b2..3d07f70 100644 --- a/KinMS.pro +++ b/KinMS.pro @@ -223,7 +223,7 @@ function kinms_samplefromarbdist_onesided,sbrad,sbprof,nsamps,seed,r_flat=r_flat theta=acos(zpos/r_3d) ;; angle out of plane xpos=((r_3d*cos(phi)*sin(theta))) ;; go to cartesian ypos=((r_3d*sin(phi)*sin(theta))) ;; go to cartesian - return,TRANSPOSE([[xpos],[ypos],[zpos]]) ;; return INCLOUDS format + return,[[xpos],[ypos],[zpos]] ;; return INCLOUDS format end function kinms_create_velfield_onesided,velrad,velprof,r_flat,inc,posang,gassigma,seed,xpos,ypos,vphasecent,vposang=vposang,vradial=vradial,inc_rad=inc_rad,posang_rad=posang_rad @@ -320,9 +320,9 @@ pro KinMS,xs,ys,vs,dx,dy,dv,beamsize,inc,gassigma=gassigma,sbprof=sbprof,sbrad=s endif ;;;; set up INCLOUDS for use - xpos=reform(INCLOUDS[0,*]/dx) - ypos=reform(INCLOUDS[1,*]/dy) - zpos=reform(INCLOUDS[2,*]/dx) + xpos=reform(INCLOUDS[*,0]/dx) + ypos=reform(INCLOUDS[*,1]/dy) + zpos=reform(INCLOUDS[*,2]/dx) if n_elements(r_flat) eq 0 then r_flat=sqrt(((xpos-(phasecen[0]/dx))^2) + ((ypos-(phasecen[1]/dy))^2)) if keyword_set(VLOS_CLOUDS) then los_vel=VLOS_CLOUDS From a873886d8b673a0028e39209ca9fa2db67680720 Mon Sep 17 00:00:00 2001 From: TimothyADavis Date: Wed, 28 Oct 2015 15:47:05 +0000 Subject: [PATCH 5/5] Updated testsuite for INCLOUDS change --- KinMS.pro | 3 +-- KinMS_testsuite.pro | 12 +++++++----- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/KinMS.pro b/KinMS.pro index 3d07f70..ceefee5 100644 --- a/KinMS.pro +++ b/KinMS.pro @@ -240,8 +240,7 @@ function kinms_create_velfield_onesided,velrad,velprof,r_flat,inc,posang,gassigm if n_elements(vradial) gt 1 then vradial_rad=interpol(vradial,velrad,r_flat) else vradial_rad=fltarr(n_elements(r_flat))+vradial if n_elements(vposang) gt 1 then vposang_rad=interpol(vposang,velrad,r_flat) else vposang_rad=vposang ;;;; - - ;;;;; include random motions + ;;;;; include random motions if keyword_set(gassigma) then begin veldisp=randomn(seed[3],n_elements(xpos)) if n_elements(gassigma) gt 1 then veldisp=temporary(veldisp)*interpol(gassigma,velrad,r_flat) else veldisp=temporary(veldisp)*gassigma diff --git a/KinMS_testsuite.pro b/KinMS_testsuite.pro index b8428e9..38b7dc6 100644 --- a/KinMS_testsuite.pro +++ b/KinMS_testsuite.pro @@ -29,10 +29,12 @@ pro KinMStest_ngc4324,_extra=_extra x1=x fx = gaussian(x,[0.1,21,3]) vel=interpol([0,50,100,175,175,175,175],[0.01,1,3,5,7,10,200],x) + phasecen=[-1,-1] + voffset=-10 ;;; ;;; Run KinMS - KinMS,xsize,ysize,vsize,dx,dy,dv,beamsize,inc,sbprof=fx,sbrad=x,velrad=x1,velprof=vel,diskthick=dickthick,nsamps=nsamps,cubeout=f,posang=posang,intflux=27.2,_extra=_extra + KinMS,xsize,ysize,vsize,dx,dy,dv,beamsize,inc,sbprof=fx,sbrad=x,velrad=x1,velprof=vel,diskthick=dickthick,nsamps=nsamps,cubeout=f,posang=posang,intflux=27.2,_extra=_extra,phasecen=phasecen,voffset=voffset ;;; ;;; Create rough moment zero and one maps and a major axis PVD from @@ -183,7 +185,7 @@ pro KinMStest_inclouds,_extra=_extra ;;;; ;;;; define where clouds are in each dimension (x,y,z) ;;;; -inclouds=[[ 40.0000, 0.00000, 0.00000],[ 39.5075, 6.25738, 0.00000],[ 38.0423, 12.3607, 0.00000],[ 35.6403, 18.1596, 0.00000],[ 32.3607, 23.5114, 0.00000],[ 28.2843, 28.2843, 0.00000],[ 23.5114, 32.3607, 0.00000],[ 18.1596, 35.6403, 0.00000],[ 12.3607, 38.0423, 0.00000],[ 6.25737, 39.5075, 0.00000],[ 0.00000, 40.0000, 0.00000],[-6.25738, 39.5075, 0.00000],[-12.3607, 38.0423, 0.00000],[-18.1596, 35.6403, 0.00000],[-23.5114, 32.3607, 0.00000],[-28.2843, 28.2843, 0.00000],[-32.3607, 23.5114, 0.00000],[-35.6403, 18.1596, 0.00000],[-38.0423, 12.3607, 0.00000],[-39.5075, 6.25738, 0.00000],[-40.0000, 0.00000, 0.00000],[-39.5075,-6.25738, 0.00000],[-38.0423,-12.3607, 0.00000],[-35.6403,-18.1596, 0.00000],[-32.3607,-23.5114, 0.00000],[-28.2843,-28.2843, 0.00000],[-23.5114,-32.3607, 0.00000],[-18.1596,-35.6403, 0.00000],[-12.3607,-38.0423, 0.00000],[-6.25738,-39.5075, 0.00000],[ 0.00000,-40.0000, 0.00000],[ 6.25738,-39.5075, 0.00000],[ 12.3607,-38.0423, 0.00000],[ 18.1596,-35.6403, 0.00000],[ 23.5114,-32.3607, 0.00000],[ 28.2843,-28.2843, 0.00000],[ 32.3607,-23.5114, 0.00000],[ 35.6403,-18.1596, 0.00000],[ 38.0423,-12.3607, 0.00000],[ 39.5075,-6.25737, 0.00000],[ 15.0000, 15.0000, 0.00000],[-15.0000, 15.0000, 0.00000],[-19.8504,-2.44189, 0.00000],[-18.0194,-8.67768, 0.00000],[-14.2856,-13.9972, 0.00000],[-9.04344,-17.8386, 0.00000],[-2.84630,-19.7964, 0.00000],[ 3.65139,-19.6639, 0.00000],[ 9.76353,-17.4549, 0.00000],[ 14.8447,-13.4028, 0.00000],[ 18.3583,-7.93546, 0.00000],[ 19.9335,-1.63019, 0.00000]] +inclouds=transpose([[ 40.0000, 0.00000, 0.00000],[ 39.5075, 6.25738, 0.00000],[ 38.0423, 12.3607, 0.00000],[ 35.6403, 18.1596, 0.00000],[ 32.3607, 23.5114, 0.00000],[ 28.2843, 28.2843, 0.00000],[ 23.5114, 32.3607, 0.00000],[ 18.1596, 35.6403, 0.00000],[ 12.3607, 38.0423, 0.00000],[ 6.25737, 39.5075, 0.00000],[ 0.00000, 40.0000, 0.00000],[-6.25738, 39.5075, 0.00000],[-12.3607, 38.0423, 0.00000],[-18.1596, 35.6403, 0.00000],[-23.5114, 32.3607, 0.00000],[-28.2843, 28.2843, 0.00000],[-32.3607, 23.5114, 0.00000],[-35.6403, 18.1596, 0.00000],[-38.0423, 12.3607, 0.00000],[-39.5075, 6.25738, 0.00000],[-40.0000, 0.00000, 0.00000],[-39.5075,-6.25738, 0.00000],[-38.0423,-12.3607, 0.00000],[-35.6403,-18.1596, 0.00000],[-32.3607,-23.5114, 0.00000],[-28.2843,-28.2843, 0.00000],[-23.5114,-32.3607, 0.00000],[-18.1596,-35.6403, 0.00000],[-12.3607,-38.0423, 0.00000],[-6.25738,-39.5075, 0.00000],[ 0.00000,-40.0000, 0.00000],[ 6.25738,-39.5075, 0.00000],[ 12.3607,-38.0423, 0.00000],[ 18.1596,-35.6403, 0.00000],[ 23.5114,-32.3607, 0.00000],[ 28.2843,-28.2843, 0.00000],[ 32.3607,-23.5114, 0.00000],[ 35.6403,-18.1596, 0.00000],[ 38.0423,-12.3607, 0.00000],[ 39.5075,-6.25737, 0.00000],[ 15.0000, 15.0000, 0.00000],[-15.0000, 15.0000, 0.00000],[-19.8504,-2.44189, 0.00000],[-18.0194,-8.67768, 0.00000],[-14.2856,-13.9972, 0.00000],[-9.04344,-17.8386, 0.00000],[-2.84630,-19.7964, 0.00000],[ 3.65139,-19.6639, 0.00000],[ 9.76353,-17.4549, 0.00000],[ 14.8447,-13.4028, 0.00000],[ 18.3583,-7.93546, 0.00000],[ 19.9335,-1.63019, 0.00000]]) ;;;; vlos_clouds=FLTARR(52) ;; in this unrealistic case each cloud is at the systemic velocity @@ -253,7 +255,7 @@ x1=a*exp(b*t)*cos(t) y1=a*exp(b*t)*sin(t) x2=a*(-1)*exp(b*t)*cos(t) y2=a*(-1)*exp(b*t)*sin(t) -inclouds=TRANSPOSE([[x1,x2],[y1,y2],[y1*0.0,y2*0.0]]) +inclouds=[[x1,x2],[y1,y2],[y1*0.0,y2*0.0]] ;;;; @@ -318,7 +320,7 @@ pro KinMStest_infits,_extra=_extra ;;;; -;;;; Read in ths FITS file and create the INCLOUDS variables based on it ;;;; +;;;; Read in the FITS file and create the INCLOUDS variables based on it ;;;; phasecent=[88,61] ;; point we wish to correspond to the phase centre in the simulation fin=readfits("test_suite/NGC1437A_FUV.fits",hdr) s=size(fin) @@ -329,7 +331,7 @@ pro KinMStest_infits,_extra=_extra index=array_indices(fin,w) x=reform(xvec[index[0,*]]) y=reform(yvec[index[1,*]]) - inclouds=TRANSPOSE([[x],[y],[y*0.0]]) + inclouds=[[x],[y],[y*0.0]] vlos_clouds=interpol([-200,0,200],[-60,0,60],y) ;; impose a flat velocity profile ;;;;