par
# compute and plot Parallactic angle
# 
# Filippenko 1982 PASP 94 715
# by GA 25/02/2003 
# type "par" to start
              # write standard  If you are at latitude=0, please move a bit...
              # define emisfero ? {Northern emisphere "N", Southern emisphere "S"}
               expand 1.5
               ctype black
               #data "parang.txt"
               data "/raid/users/altavilla/public_html/tmp/parang.txt"
               device postencap "/raid/users/altavilla/public_html/tmp/prova.ps"

		#define long ? {observer longitudine ( + West, - Est; decimal degree)}
	        #define fi ? {observer latitude (decimal degree, only + values!)(+45.83 Ekar, +42.26 Loiano, +33.36 Palomar, +25.76 TNG, +37.22 CAHA, -29.25 La Silla, -24.67 Paranal, -29.0 Las Campanas)}


               read {dati 1}
              
   
              define fi1 (dati[0])
              define fi2 (dati[1])
              define fi3 (dati[2])

		define fi ($fi1+($fi2+$fi3/60)/60)
   
	       if($fi>0){parnord}
	       if($fi<0){parsud}	



parnord

		

               define pi (acos(-1))
		
                define RA1 (dati[6])
                define RA2 (dati[7])
                define RA3 (dati[8])

                define RA ($RA1+($RA2+$RA3/60)/60) #RA decimal
                write standard RA=$RA            
 
                define ST1  (dati[3])
                define ST2  (dati[4])
                define ST3  (dati[5])

                define ST ($ST1+($ST2+$ST3/60)/60) #ST decimal
                write standard ST=$ST



               
                define h0 ($ST-$RA)
                if($h0<-12){define h0 ($h0+24)}
                if($h0>+12){define h0 ($h0-24)}

 
		#define h0 ? {hour angle ? (decimal hours + west, - est)}
		#define long ? {observer longitudine ( + West, - Est; decimal degree) (-0.76 Ekar)}
               



 

                define delta1 (dati[9])#target declination (this value must be different from the latitude!!!)
                define delta2 (dati[10])
                define delta3 (dati[11])

                define delta ($delta1+($delta2+$delta3/60)/60) #delta decimal
                write standard delta=$delta


                
                define h (abs($h0))
		define hrad ($h*15/(180/$pi))
		define firad ($fi/(180/$pi))
                define deltarad ($delta/(180/$pi))
		define deltamin (dati[12]) #{Minimum declination value to be plotted? (-40)}
      		define deltamax (dati[13]) #{Maximum declination value to be plotted? (<90)}
     	        define step (dati[14]) #{Declination step? (10)}        

                
                erase
                limits 0 8 0 180
                ctype black
                box
                xlabel H.A.
                ylabel P.A.
                
                set xline=0,8,.1
                set yline1=xline*0
                set yline2=xline*0
                set yline3=xline*0
                define n 0
	        do hh=0,8,.1 {
		
                	
			define hrad ($hh*15/(180/$pi))
			define firad ($fi/(180/$pi))
              	        define deltarad ($delta/(180/$pi))
		
			define etarad (asin((sin($hrad)*cos($firad))/sqrt(1-(sin($firad)*sin($deltarad)+cos($firad)*cos($deltarad)*cos($hrad))**2)))

			define eta1 ($etarad*(180/$pi))
            	        define eta2 (180-$eta1)
		 	define eta3 (-180-$eta1)
			

                        set x=$hh
                        
                        set yline1[$n]=$eta1
                        set yline2[$n]=$eta2
                        set yline3[$n]=$eta3
                            
                       define n ($n+1)
	                }
                        
                ctype blue
                lweight 3
                if($delta<$fi){
                               set yline4= yline1
                               connect xline yline4
                               }
                
                if($delta>$fi){
                              set xlined=xline
                              set diff= yline2-yline1   
                              sort {diff xlined}
                              define  x1 (xlined[0])
                              set yline4= xline<$x1 ? yline2 : yline1 
                              connect xline  yline4
                               }
                if($delta==$fi){write standard *******************************************************
                               write standard CAUTION!!!!
                               write standard if Dec=$delta=Latitude=$fi
                               write standard ambigous angle
                               write standard *******************************************************
                               relocate (7000 20000)
                               ctype red
                               label Please enter declination 
                               relocate (7000 19000) 
                               label different from the latitude!
                               ctype black
                               }

              interp xline yline4 $h PARALLACTIC_ANGLE
        
              if($h0<0) {set PARALLACTIC_ANGLE=- PARALLACTIC_ANGLE}     

		write standard ************************************************
                write standard
		write standard latitude=$fi (decimal degree)
		write standard dec=$delta   (decimal degree)
                 write standard HOUR ANGLE=$h0
		set Hour_Angle = ($h0)
                               print {PARALLACTIC_ANGLE}
                 write standard WARNING: AFOSC slit is oriented  E-W (and not N-S as usual) so: 
                 #if(PARALLACTIC_ANGLE>180) {set AFOSC_PA=PARALLACTIC_ANGLE-180-90}  180<P.A<180
                 if(PARALLACTIC_ANGLE<180&&PARALLACTIC_ANGLE>0) {set AFOSC_PA=PARALLACTIC_ANGLE-90}  
                 if(PARALLACTIC_ANGLE<0&&PARALLACTIC_ANGLE>-180) {set AFOSC_PA=PARALLACTIC_ANGLE+180-90}  
                 #if(PARALLACTIC_ANGLE<-180){ set AFOSC_PA=PARALLACTIC_ANGLE-180-90}                     
                 print {AFOSC_PA}
                 print "/raid/users/altavilla/public_html/tmp/PA.file" {Hour_Angle}
                 print +"/raid/users/altavilla/public_html/tmp/PA.file" {PARALLACTIC_ANGLE}
		 print +"/raid/users/altavilla/public_html/tmp/PA.file" {AFOSC_PA} 
		write standard ************************************************



                  set xlab=4
                  set delta22=$delta
                  set ynome=string(delta22)
                  ptype ynome
                  interp xline yline4 xlab  ylab
                  set ylab=ylab+5
                  points xlab ylab     


########################################### plot for different values of declination  ###################################

       lweight 1  

       do  delta2=$deltamin,$deltamax,$step { 
                if($delta2==$fi){define delta2 ($fi+0.00001)
                                write standard =======================================
                                write standard ATTENTION:
                                write standard if  declination = latitude
                                write standard not reliable values!!!!!!!!!!!!!!!!!!
                                write standard =======================================}

                define delta $delta2         
                set xline=0,8,.1
                set yline1=xline*0
                set yline2=xline*0
                set yline3=xline*0
                
                define n 0
	        do hh=0,8,.1 {
                	define hh (abs($hh))
			define hrad ($hh*15/(180/$pi))
			define firad ($fi/(180/$pi))
              	        define deltarad ($delta/(180/$pi))
		
			define etarad (asin((sin($hrad)*cos($firad))/sqrt(1-(sin($firad)*sin($deltarad)+cos($firad)*cos($deltarad)*cos($hrad))**2)))

			define eta1 ($etarad*(180/$pi))
            	        define eta2 (180-$eta1)
		 	define eta3 (-180-$eta1)
			

                        set x=$hh
                                                
                        set yline1[$n]=$eta1
                        set yline2[$n]=$eta2
                        set yline3[$n]=$eta3
                            
                        
                       define n ($n+1)
	                }


		

                ctype green
                if($delta<$fi){
                               set yline4= yline1
                               connect xline yline4
                               }
                
                if($delta>$fi){
                              set xlined=xline
                              set diff= yline2-yline1   
                              sort {diff xlined}
                              define  x1 (xlined[0])
                              set yline4= xline<$x1 ? yline2 : yline1 
                              connect xline  yline4
                               }
                if($delta==$fi){write standard *******************************************************
                               write standard CAUTION!!!!
                               write standard if Dec=$delta=Latitude=$fi
                               write standard ambigous angle
                               write standard *******************************************************
                               relocate (7000 20000)
                               ctype red
                               label Please enter declination 
                               relocate (7000 19000) 
                               label different from the latitude!
                               ctype black
                               }
                  
                  set xlab=2+(abs($delta2)/100)
                  set delta22=$delta2
                  set ynome=string(delta22)
                  ptype ynome
                  interp xline yline4 xlab  ylab
                  set ylab=ylab+5

                  points xlab ylab     

                         }       

                  write standard DONE

##################cambio emisfero

parsud

               define pi (acos(-1))
			
	        define RA1 (dati[6])# {RA hh}
                define RA2 (dati[7])# {RA mm}
                define RA3 (dati[8])# {RA ss}

                define RA ($RA1+($RA2+$RA3/60)/60) #RA decimal
                write standard RA=$RA            
 
                define ST1 (dati[3])# {ST hh}
                define ST2 (dati[4])# {ST mm}
                define ST3 (dati[5])# {ST ss}

                define ST ($ST1+($ST2+$ST3/60)/60) #ST decimal
                write standard ST=$ST



               
                define h0 ($ST-$RA)
                if($h0<-12){define h0 ($h0+24)}
                if($h0>+12){define h0 ($h0-24)}

                
		

                

                define delta1 (dati[9]) #target declination (this value must be different from the latitude!!!)
                define delta2 (dati[10])
                define delta3 (dati[11])

                define delta ($delta1+($delta2+$delta3/60)/60) #delta decimal
                write standard delta=$delta




                define h (abs($h0))
		define hrad ($h*15/(180/$pi))
		define firad ($fi/(180/$pi))
                define deltarad ($delta/(180/$pi))
                define deltamin (dati[12]) #{Minimum declination value to be plotted? (>-90)}
      		define deltamax (dati[13])#{Maximum declination value to be plotted? (+40)}
     	        define step (dati[14]) #{Declination step? (10)}        

              
                
                erase
                limits 0 8 0 180
                ctype black
                box
                xlabel A.H.
                ylabel P.A.
                
                set xline=0,8,.1
                set yline1=xline*0
                set yline2=xline*0
                set yline3=xline*0
                define n 0
	        do hh=0,8,.1 {
		
                	
			define hrad ($hh*15/(180/$pi))
			define firad ($fi/(180/$pi))
              	        define deltarad ($delta/(180/$pi))
		
			define etarad (asin((sin($hrad)*cos($firad))/sqrt(1-(sin($firad)*sin($deltarad)+cos($firad)*cos($deltarad)*cos($hrad))**2)))

			define eta1 ($etarad*(180/$pi))
            	        define eta2 (180-$eta1)
		 	define eta3 (-180-$eta1)
			

                        set x=$hh
                        
                        set yline1[$n]=$eta1
                        set yline2[$n]=$eta2
                        set yline3[$n]=$eta3
                            
                       define n ($n+1)
	                }
                       #ptype 3 0 
                       #ctype red
                       #points xline yline1
                       #ctype green
                       #points xline yline2
                       #ctype cyan
                       #points xline yline3

                ctype blue
                lweight 3
                if($delta>$fi){ 
                               set yline4= yline2
                               connect xline yline4
                               }
                
                if($delta<$fi){
                              set xlined=xline
                              set diff= yline2-yline1   
                              sort {diff xlined}
                              define  x1 (xlined[0])
                              set yline4= xline<$x1 ? yline1 : yline2 
                              connect xline  yline4
                               }
                if($delta==$fi){write standard *******************************************************
                               write standard CAUTION!!!!
                               write standard if Dec=$delta=Latitude=$fi
                               write standard ambigous angle
                               write standard *******************************************************
                               relocate (7000 20000)
                               ctype red
                               label Please enter declination 
                               relocate (7000 19000)
                               label different from the latitude!
                               ctype black
                               }

              interp xline yline4 $h PARALLACTIC_ANGLE
        
              if($h0<0) {set PARALLACTIC_ANGLE=- PARALLACTIC_ANGLE}     

		write standard ************************************************
                write standard
		write standard latitude=$fi (decimal degree)
		write standard dec=$delta   (decimal degree)
                write standard HOUR ANGLE=$h0
                set Hour_Angle = ($h0)
                 print {PARALLACTIC_ANGLE}
                 print "/raid/users/altavilla/public_html/tmp/PA.file" {Hour_Angle}
                 print +"/raid/users/altavilla/public_html/tmp/PA.file" {PARALLACTIC_ANGLE}
		
		write standard ************************************************



                  set xlab=4
                  set delta22=$delta
                  set ynome=string(delta22)
                  ptype ynome
                  interp xline yline4 xlab  ylab
                  set ylab=ylab+5
                  points xlab ylab     



########################################### general plot ###################################

       lweight 1  

       do  delta2=$deltamin,$deltamax,$step { 
                if($delta2==$fi){define delta2 ($fi+0.00001)
                                write standard =======================================
                                write standard ATTENTION:
                                write standard if  declination = latitude
                                write standard not reliable values!!!!!!!!!!!!!!!!!!
                                write standard =======================================}

                define delta $delta2         
                set xline=0,8,.1
                set yline1=xline*0
                set yline2=xline*0
                set yline3=xline*0
                
                define n 0
	        do hh=0,8,.1 {
                	define hh (abs($hh))
			define hrad ($hh*15/(180/$pi))
			define firad ($fi/(180/$pi))
              	        define deltarad ($delta/(180/$pi))
		
			define etarad (asin((sin($hrad)*cos($firad))/sqrt(1-(sin($firad)*sin($deltarad)+cos($firad)*cos($deltarad)*cos($hrad))**2)))

			define eta1 ($etarad*(180/$pi))
            	        define eta2 (180-$eta1)
		 	define eta3 (-180-$eta1)
			

                        set x=$hh
                                                
                        set yline1[$n]=$eta1
                        set yline2[$n]=$eta2
                        set yline3[$n]=$eta3
                            
                        
                       define n ($n+1)
	                }


		

                ctype green
                if($delta>$fi){
                               set yline4= yline2
                               connect xline yline4
                               }
                
                if($delta<$fi){
                              set xlined=xline
                              set diff= yline2-yline1   
                              sort {diff xlined}
                              define  x1 (xlined[0])
                              set yline4= xline<$x1 ? yline1 : yline2 
                              connect xline  yline4
                               }
                if($delta==$fi){write standard *******************************************************
                               write standard CAUTION!!!!
                               write standard if Dec=$delta=Latitude=$fi
                               write standard ambigous angle
                               write standard *******************************************************
                               relocate (7000 20000)
                               ctype red
                               label Please enter declination 
                               relocate (7000 19000)
                               label different from the latitude!
                               ctype black
                               }
                  
                  set xlab=2+(abs($delta2)/100)
                  set delta22=$delta2
                  set ynome=string(delta22)
                  ptype ynome
                  interp xline yline4 xlab  ylab
                  set ylab=ylab+5

                  points xlab ylab     

                         }       
                  write standard DONE
